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Abstract 


The tools of weakly coupled phase oscillator theory have had a profound impact on the 
neuroscience community, providing insight into a variety of network behaviours ranging from 
central pattern generation to synchronisation, as well as predicting novel network states such 
as chimeras. However, there are many instances when this theory is expected to break down, 
say in the presence of strong coupling, or must be carefully interpreted, as in the presence of 
stochastic forcing. There are also surprises in the dynamical complexity of the attractors that 
can robustly appear - for example, heteroclinic network attractors. In this review we present 
a set of mathematical tools that are suitable for addressing the dynamics of oscillatory neural 
networks, broadening from a standard phase oscillator perspective to provide a practical frame¬ 
work for further successful applications of mathematics to understanding network dynamics in 
neuroscience. 
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1 Introduction 

Coupled oscillator theory is now a pervasive part of the theoretical neuroscientist’s toolkit for study¬ 
ing the dynamics of models of biological neural networks. Undoubtedly this technique originally 
arose in the broader scientific community through a fascination with understanding synchronisa¬ 
tion in networks of interacting heterogeneous oscillators, and can be traced back to the work of 
Huygens on “an odd kind of sympathy” between coupled pendulum clocks” |143j . Subsequently the 
theory has been developed and applied to the interaction between organ pipes |268j , phase-locking 
phenomena in electronic circuits [257], the analysis of brain rhythms [346] . chemical oscillations 
|I73|, cardiac pacemakers |218j . circadian rhythms |199j . flashing fireflies |98j . coupled Josephson 
junctions [347] . rhythmic applause [230] . animal flocking [131] . fish schooling [246] . and behaviours 
in social networks |33|. For a recent overview of the application of coupled phase oscillator theory 
to areas as diverse as vehicle coordination, electric power networks, and clock synchronisation in 
decentralised networks see the recent survey article by Dorfler and Bullo |93j . 

Given the widespread nature of oscillations in neural systems it is no surprise that the science 
of oscillators has found such ready application in neuroscience |334j . This has proven especially 
fruitful for shedding light on the functional role that oscillations can play in feature binding |2Ml 
1329] , cognition [341] , memory processes |110j , odour perception |345l 1187] , information transfer 
mechanisms m, inter-limb coordination UMlCSTj, and the generation of rhythmic motor output 
IMI. Neural oscillations also play an important role in many neurological disorders, such as 
excessive synchronisation during seizure activity in epilepsy |220L I76j . tremor in patients with 
Parkinson’s disease |324] or disruption of cortical phase synchronisation in schizophrenia |l6]. As 
such it has proven highly beneficial to develop methods for the control of (de)synchronisation in 
oscillatory networks, as exemplified by the work of Tass et al. [315] 1316] for therapeutic brain 
stimulation techniques. From a transformative technological perspective, oscillatory activity is 
increasingly being used to control external devices in brain-computer interfaces, in which subjects 
can control an external device by changing the amplitude of a particular brain rhythm |267] . 

Neural oscillations can emerge in a variety of ways, including intrinsic mechanisms within in¬ 
dividual neurons or by interactions between neurons. At the single neuron level, sub-threshold 
oscillations can be seen in membrane voltage as well as rhythmic patterns of action potentials. 
Both can be modelled using the Hodgkin-Huxley conductance formalism, and analysed mathemat¬ 
ically with dynamical systems techniques to shed light on the mechanisms that underly various 
forms of rhythmic behaviour, including tonic spiking and bursting (see e.g. |73|). The high dimen¬ 
sionality of biophysically realistic single neuron models has also encouraged the use of reduction 
techniques, such as the separation of time-scales recently reviewed in [272L 1172] . or the use of 
phenomenological models, such as FitzHugh-Nagumo (FHN) |113] . to regain some level of mathe¬ 
matical tractability. This has proven especially useful when studying the response of single neurons 
to forcing itself a precursor to understanding how networks of interacting neurons can behave. 
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When mediated by synaptic interactions, the repetitive firing of pre-synaptic neurons can cause 
oscillatory activation of post-synaptic neurons. At the level of neural ensembles, synchronised ac¬ 
tivity of large numbers of neurons gives rise to macroscopic oscillations, which can be recorded with 
a micro-electrode embedded within neuronal tissue as a voltage change referred to as a local field 
potential (LFP). These oscillations were first observed outside the brain by Hans Berger in 1924 
m in electroencephalogram (EEG) recordings, and have given rise to the modern classification of 
brain rhythms into frequency bands for alpha activity (8-13 Hz) (recorded from the occipital lobe 
during relaxed wakefulness), delta (1-4 Hz), theta (4-8 Hz), beta (13-30 Hz) and gamma (30-70 Hz). 
The latter rhythm is often associated with cognitive processing, and it is now common to link large 
scale neural oscillations with cognitive states, such as awareness and consciousness. For example, 
from a practical perspective the monitoring of brain states via EEG is used to determine depth of 
anaesthesia |336j . Such macroscopic signals can also arise from interactions between different brain 
areas, the thalamo-cortical loop being a classic example |286j . Neural mass models (describing 
the coarse grained activity of large populations of neurons and synapses) have proven especially 
useful in understanding EEG rhythms [356], as well as in augmenting the dynamic causal modelling 
framework (driven by large scale neuroimaging data) for understanding how event-related responses 
result from the dynamics of coupled neural populations |H3| . 

One very influential mathematical technique for analysing networks of neural oscillators, whether 
they be built from single neuron or neural mass models, has been that of weakly coupled oscilla¬ 
tor theory, as comprehensively described by Hoppensteadt and Izhikevich |142) . In the limit of 
weak coupling between limit cycle oscillators invariant manifold theory m and averaging theory 
|128| can be used to reduce the dynamics to a set of phase equations in which the relative phase 
between oscillators is the relevant dynamical variable. This approach has been applied to neural 
behaviour ranging from that seen in small rhythmic networks |161) up to the whole brain [60] . 
Despite the powerful tools and wide-spread use afforded by this formalism, it does have a number 
of limitations (such as assuming the persistence of the limit cycle under coupling) and it is well 
to remember that there are other tools from the mathematical sciences relevant to understanding 
network behaviour. In this review we wrap the weakly coupled oscillator formalism in a variety of 
other techniques ranging from symmetric bifurcation theory and groupoid formalisms through to 
more “physics-based” approaches for obtaining reduced models of large networks. This highlights 
the regimes where the standard formalism is applicable, and provides a set of complementary tools 
when it does not. These are especially useful when investigating systems with strong coupling, or 
ones for which the rate of attraction to a limit cycle is slow. 

In § [^ we review some of the key mathematical models of oscillators in neuroscience, ranging 
from single neuron to neural mass, as well as introduce the standard machinery for describing 
synaptic and gap junction coupling. We then present in § [^ an overview of some of the more 
powerful mathematical approaches to understanding the collective behaviour in coupled oscillator 
networks, mainly drawn from the theory of symmetric dynamics. We touch upon the master 
stability function approach and the groupoid formalism for handling coupled cell systems. In § [^ 
we review some special cases where it is either possible to say something about the stability of the 
globally synchronous state in a general setting, or that of phase-locked states for strongly coupled 
networks of integrate-and-fire neurons. The challenge of the general case is laid out in § [^ where 
we advocate the use of phase-amplitude coordinates as a starting point for either direct network 
analysis or network reduction. To highlight the importance of dynamics off cycle we discuss the 
phenomenon of shear-induced chaos. In the same section we review the reduction to the standard 
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phase-only description of an oscillator, covering the well known notions of isochrons and phase 
response curves. The construction of phase interaction functions for weakly coupled phase oscillator 
networks is covered in § [^ together with tools for analysing phase-locked states. Moreover, we go 
beyond standard approaches and describe the emergence of turbulent states in continuum models 
with non-local coupling. Another example of something more complicated than a periodic attractor 
is that of a heteroclinic attractor, and these are the subject of §[^ The subtleties of phase reduction 
in the presence of stochastic forcing are outlined in § [^ The search for reduced descriptions of very 
large networks is the topic of § |^ where we cover recent results for Winfree networks that provide 
an exact mean-field description in terms of a complex order parameter. This approach makes use 
of the Ott-Antonsen ansatz that has also found application to chimera states, and which we discuss 
in a neural context. In §[^we briefly review some examples where the mathematics of this review 
have been applied, and finally in § we discuss some of the many open challenges in the field of 
neural network dynamics. 

We will assume the reader has familiarity with the following: 

• The basics of nonlinear differential equation descriptions of dynamical systems such as linear 
stability and phase-plane analysis. 

• Ideas from the qualitative theory of differential equations/dynamical systems such as asymp¬ 
totic stability, attractors and limit cycles. 

• Generic codimension-one bifurcation theory of equilibria (saddle-node, Hopf) and of periodic 
orbits (saddle-node of limit cycles, heteroclinic, torus, flip). 

There are a number of texts that cover this material very well in the context of neuroscience 
modelling, for example |146l I105j . At the end we include a glossary of some abbreviations that are 
introduced in the text. 

2 Neurons and neural populations as oscillators 

Nonlinear ionic currents, mediated by voltage-gated ion channels, play a key role in generating 
membrane potential oscillations and action potentials. There are many ordinary differential equa¬ 
tion (ODE) models for voltage oscillations, ranging from biophysically detailed conductance-based 
models through to simple integrate-and-fire (IF) caricatures. This style of modelling has also proved 
fruitful at the population level, for tracking the averaged activity of a near synchronous state. In 
all these cases bifurcation analysis is especially useful for classifying the types of oscillatory (and 
possibly resonant) behaviour that are possible. Here we give a brief overview of some of the key 
oscillator models encountered in computational neuroscience, as well as models for electrical and 
chemical coupling necessary to build networks. 

2.1 The Hodgkin-Huxley model and its planar reduction 

The work of Hodgkin and Huxley in elucidating the mechanism of action potentials in the squid 
giant axon is one of the major breakthroughs of dynamical modelling in physiology [137] . and see 
m for a review. Their work underpins all modern electrophysiological models, exploiting the 
observation that cell membranes behave much like electrical circuits. The basic circuit elements are 
1) the phospholipid bilayer, which is analogous to a capacitor in that it accumulates ionic charge as 
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the electrical potential across the membrane changes; 2) the ionic permeabilities of the membrane, 
which are analogous to resistors in an electronic circuit; and 3) the electrochemical driving forces, 
which are analogous to batteries driving the ionic currents. These ionic currents are arranged in 
a parallel circuit. Thus the electrical behaviour of cells is based upon the transfer and storage of 
ions such as K"*" and Na"*". 

Our goal here is to illustrate, by exploiting specific models of excitable membrane, some of the 
concepts and techniques which can be used to understand, predict, and interpret the excitable and 
oscillatory behaviours that are commonly observed in single cell electrophysiological recordings. 
We begin with the mathematical description of the Hodgkin-Huxley model. 

The standard dynamical system for describing a neuron as a spatially isopotential cell with 
constant membrane potential V is based upon conservation of electric charge, so that 

C^V = lion + I, 
at 

where C is the cell capacitance, I the applied current and lion represents the sum of individual 
ionic currents: 

/ion = -gxiv - Vk) - gNa{V - Vno) " g^V - Vl). 

In the Hodgkin-Huxley (HH) model the membrane current arises mainly through the conduction of 
sodium and potassium ions through voltage dependent channels in the membrane. The contribution 
from other ionic currents is assumed to obey Ohm’s law (and is called the leak current). The gx, 
gxa and gi are conductances (conductance=l/resistance). 

The great insight of Hodgkin and Huxley was to realise that gx depends upon four activation 
gates: gx = gx'^^'^i whereas gxa depends upon three activation gates and one inactivation gate: 
gNa = gNa'^^^h. Here the gating variables all obey equations of the form 


d X^{V)-m 

dt rx{V) ’ 


X G {m, n, h}. 


The conductance variables described by X take values between 0 and 1 and approach the asymptotic 
values Xoo{V) with time constants rx(H). These six functions are obtained from fits with experi¬ 
mental data. It is common practice to write tx{V) = 1/{ax{y) + f^xiV))-, X^oiV) = ax{V)Tx{y), 
where a and /3 have the interpretation of opening and closing channel transition rates respectively. 
The details of the HH model are provided in Appendix A for completeness. A numerical bifurcation 
diagram of the model in response to constant current injection is shown in Fig. illustrating that 
oscillations can emerge in a Hopf bifurcation with increasing drive. 

The mathematical forms chosen by Hodgkin and Huxley for the functions tx and X^o, X G 
{m, n, h}, are all transcendental functions. Both this and the high dimensionality of the model make 
analysis difficult. However, considerable simplification is attained with the following observations: 
(i) Tm{V) is small for all V so that the variable m rapidly approaches its equilibrium value moo{V), 
and (ii) the equations for h and n have similar time-courses and may be slaved together. This has 
been put on a more formal footing by Abbott and Kepler [1], to obtain a reduced planar model 
for (V, U) obtained from the full Hodgkin-Huxley model under the replacement m —)> moo(K) and 
X = Xoo{U) for X G {n, h} with a prescribed choice of U. The phase-plane and nullclines for this 
system are shown in Fig. 

For zero external input the fixed point is stable and the neuron is said to be excitable. When a 
positive external current is applied the low-voltage portion of the V nullcline moves up whilst the 
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Figure 1: Bifurcation diagram of the Hodgkin-Huxley model as a function of the external drive I. 
The green circles show the amplitude of a stable limit cycle and the blue circles indicate unstable 
limit cycle behaviour. The solid red line shows the stable fixed point and the black line shows the 
unstable fixed point. Details of the model are in Appendix A. 



Figure 2: Nullclines (red for V and green for U) of the reduced HH neuron mode, obtained using 
the reduction technique of Abbott and Kepler [T], for the oscillatory regime (/ = 10) capable of 
supporting a periodic train of spikes. The periodic orbit is shown in blue. 
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high-voltage part remains relatively unchanged. For sufficiently large constant external input the 
intersection of the two nullclines falls within the portion of the V nullcline with positive slope. In 
this case the fixed point is unstable and the system may support a limit cycle. The system is said to 
be oscillatory as it may produce a train of action potentials. Action potentials may also be induced 
in the absence of an external current for synaptic stimuli of sufficient strength and duration. This 
simple planar model captures all of the essential features of the original HH model yet is much 
easier to understand from a geometric perspective. Indeed the model is highly reminiscent of the 
famous FHN model, in which the voltage nullcline is taken to be a cubic function. Both models 
show the onset of repetitive firing at a non-zero frequency as observed in the HH model (when 
an excitable state loses stability via a subcritical Hopf bifurcation). However, unlike real cortical 
neurons they cannot fire at arbitrarily low frequency. This brings us to consider modifications of 
the original HH formalism to accommodate bifurcation mechanisms from excitable to oscillatory 
behaviours that can respect this experimental observation. 

2.2 The cortical model of Wilson 

Many of the properties of real cortical neurons can be captured by making the equation for the 
recovery variable of the FHN equations quadratic (instead of linear). We are thus led to the cortical 
model of Wilson |348j : 

C-^v = f(v) - w + I, f(v) = v(a-v)(v-l), 

—w = /3(v — vi)(v — V 2 ) — 'jw, 

where 0 < a < 1 and 6, 7 , /r > 0. Here v is like the membrane potential V, and w plays the role of a 
gating variable. In addition to the single fixed point of the FHN model it is possible to have another 
pair of fixed points, as shown in Fig. (right). As I increases two fixed points can annihilate at 
a saddle node on an invariant curve (SNIC) bifurcation at / = Ic |146| . In the neighbourhood of 
this global bifurcation the firing frequency scales like \/l — Ic- For large enough I there is only one 
fixed point on the middle branch of the cubic, as illustrated in Fig. (left). In this instance an 
oscillatory solution occurs via the same mechanism as for the FHN model. 

2.3 Morris-Lecar with homoclinic bifurcation 

A SNIC bifurcation is not the only way to achieve a low firing rate in a single neurone model. 
It is also possible to achieve this via a homoclinic bifurcation, as is possible in the Morris-Lecar 
(ML) model |225| . This was originally developed as a model for barnacle muscle fiber. Morris 
and Lecar introduced a set of coupled ordinary differential equations (ODEs) incorporating two 
ionic currents: an outward going, non-inactivating potassium current and an inward going, non¬ 
inactivating calcium current. Assuming that the calcium currents operate on a much faster time 
scale than the potassium current one they formulated the following two dimensional system: 

C^V = gL{VL -V)+ gKw{VK - H) + gcamoo{V){Vca - H) + I, 









Figure 3: Phase portrait for cortical neuron model with quadratic recovery variable, a = 0.1, 
/3 = 7 = 0.5, vi = 0,V2 = 0.2. The voltage nullcline is shown in red and that of the recovery 
variable in green. Left: 1 = 0, showing a stable fixed point (black filled circle), a saddle (grey filled 
circle) and an unstable fixed point (white filled circle). Right: I = 0.1, where there is an unstable 
fixed point (white filled circle) with a stable limit cycle (in blue) for C = 0.01. 


with moo{V) = 0.5(1 + tanh[(l/ — Vi)/V 2 ])) 'U^oo(R) = 0.5(1 + tanh[(F — Vs)/!^]), and A(F) = 
(/>cosh[(F — 1 / 3 ) 7 ( 214 )]. Here w represents the fraction of K"*" channels open, and the Ca^'*' channels 
respond to V so rapidly that we assume instantaneous activation. Here gi is the leakage conduc¬ 
tance, gK,gCa are the potassium and calcium conductances, Vl, Vk, Vca are corresponding reversal 
potentials, moo(R), rfoo(F) are voltage-dependent gating functions and A(F) is a voltage-dependent 
rate. Referring to Fig. as I decreases the periodic orbit grows in amplitude, it comes closer to 
a saddle point and slows down such that near the homoclinic bifurcation, where the orbit collides 
with the saddle at I = Ic, the frequency of oscillation scales as —l/log(/ — Ic)- 


2.4 Integrate-and-fire 

Although conductance-based models like that of Hodgkin and Huxley provide a level of detail that 
helps us to understand how the kinetics of channels (with averaged activation and inactivation 
variables) can underlie action-potential generation their high dimensionality is a barrier to studies 
at the network level. The goal of a network-level analysis is to predict emergent computational 
properties in populations and recurrent networks of neurons from the properties of their component 
cells. Thus simpler (lower dimensional and hopefully mathematically tractable) models are more 
appealing - especially if they fit single neuron data. 

A one-dimensional nonlinear IF model takes the form 


= f(v) + ^ 


ifv< vth, 


( 1 ) 


such that v(t) is reset to ur just after reaching the value uth- In other words we seek a piecewise 
discontinuous solution v(t) of 0 such that 


lim u(to) = Vt\i implies that lim v{t) = ur. 

t ^ 


Firing times are defined iteratively according to 


Tn = inf{t \v{t) > uth ; t> Tn-l}. 
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Figure 4: Phase portrait of the Morris-Lecar model at / = 0.075 with C = 1, 14 = —0.7, Vl = —0.5, 
^Ca = l,gK = 2, gL = 0.5, Pi = -0.01, P2 = 0.15, gca = 1-33, P 3 = 0.1, P 4 = 0.145 and </) = 1/3. 
The voltage nullcline is shown in red and that of the gating variable in green. The hlled black circle 
indicates a stable fixed point, the grey hlled circle a saddle and the hlled white circle an unstable 
hxed point. The periodic orbit is shown in blue. 


Real cortical data can be very accurately ht using: 

f{v) = ul — u + 

with ul = —68.5 mV, r = 3.3 ms, = —61.5 mV and k = 4 mV |36]. There are many varieties 
of nonlinear IF model, with the quadratic one |185| being well known as a precursor for the planar 
Izhikevich spiking model |148j . itself capable of generating a wide variety of bring patterns, including 
bursting and chattering as well as regular spiking. For a more thorough discussion of IF models 
and the challenges of analysing non-smooth dynamical systems we refer the reader to 1771. 

2.5 Neuronal coupling 

At a synapse presynaptic bring results in the release of neurotransmitters that cause a change in 
the membrane conductance of the postsynaptic neuron. This postsynaptic current may be written 
Is{t) = gss{t){Vs — V{t)) where V{t) is the voltage of the postsynaptic neuron, Vs is the membrane 
reversal potential and gs is a constant conductance. The variable s corresponds to the probability 
that a synaptic receptor channel is in an open conducting state. This probability depends on the 
presence and concentration of neurotransmitter released by the presynaptic neuron. The sign of I 4 
relative to the resting potential (which without loss of generality we may set to zero) determines 
whether the synapse is excitatory (I 4 > 0) or inhibitory (I 4 < 0). 

The effect of some synapses can be described with a function that hts the shape of the postsy¬ 
naptic response due to the arrival of action potential at the pre-synaptic release site. A postsynaptic 
potential (PSP) s{t) would then be given by s{t) = g{t — T), t > T where T is the arrival time of 
a pre-synaptic action potential and r]{t) hts the shape of a realistic PSP (with g{t) = 0 for t < 0). 
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A common (normalised) choice for r]{t) is a difference of exponentials: 

or the alpha function obtained from ([^ in the limit (5 ^ a. The PSP arising from a train 

of action potentials is given by 

s{i) ='^v{t-Tm), (3) 

mGZ 

where denotes the arrival time of the mth action potential at a synapse. 

Interestingly even purely inhibitory synaptic interactions between non-oscillatory neurons can 
create oscillations at the network level, and can give rise to central pattern generators of “half¬ 
center” type [56]. To see this we need only consider a pair of (reduced) Hodgkin-Huxley neurons 
with mutual reciprocal inhibition mediated by an a-function synapse with a negative reversal 
potential. The phenomenon of anode break excitation (whereby a neuron fires an action potential 
in response to termination of a hyperpolarising current) can underlie a natural anti-phase rhythm, 
and is best understood in terms of the phase plane shown in Fig. In this case inhibition will 
effectively move the voltage nullcline down, and the system will equilibrate to a new hyperpolarised 
state. Upon release of inhibition the fixed point will move to a higher value, though to reach this 
new state the trajectory must jump to the right hand branch (since now dU /dt > 0). An example 
is shown in Fig. 



Figure 5: A Half Center oscillator built from two (reduced) Hodgkin-Huxley neurons with mutual 
reciprocal inhibition modelled by an a-function synapse and a negative reversal potential. The 
inset shows the voltage traces (solid and dashed lines) for the two neurons. The solid blue line in 
the {V,U) space shows the common orbit of the two neurons (though each neuron travels half a 
period out of phase with the other). Parameters: a = 1 ms, 14 = —100 mV, Qs = 100 mS cm“^. 
Spike times determined by a voltage threshold at U = —35 mV. 


Gap junctions differ from chemical synapses in that they allow for direct communication between 
cells. They are typically formed from the juxtaposition of two hemichannels (connexin proteins) 
and allow the free movement of ions or molecules across the intercellular space separating the 
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plasma membrane of one cell from another. As well as being found in the neocortex, they occur in 
many other brain regions, including the hippocampus, inferior olivary nucleus in the brain stem, 
the spinal cord, and the thalamus m- Without the need for receptors to recognise chemical 
messengers, gap junctions are much faster than chemical synapses at relaying signals. The synaptic 
delay for a chemical synapse is typically in the range 1-lOOms, while the synaptic delay for an 
electrical synapse may be only about 0.2ms. 

It is common to view the gap-junction as nothing more than a channel that conducts current 
according to a simple ohmic model. For two neurons with voltages Vi and vj the current flowing 
into cell i from cell j is given by Iga.p{vi,Vj) = g{vj — Vi), where g is the constant strength of the 
gap junction conductance. They are believed to promote synchrony between oscillators (e.g. see 
|325| ), though the story is more subtle than this as we shall discuss in § [^ 

2.6 Neural mass models 

As well as supporting oscillations at the single neuron level, brain tissue can also generate oscilla¬ 
tions at the tissue level. Rather than model this using networks built from single neuron models, it 
is has proven especially useful to develop low dimensional models to mimic the collection of thou¬ 
sands of near identical interconnected neurons with a preference to operate in synchrony. These 
are often referred to as neural mass models, with state variables that track coarse grained measures 
of the average membrane potential, firing rates or synaptic activity. They have proven especially 
useful in the description of human EEG power spectra |193j , as well as resting brain state activity 
[51] and mesoscopic brain oscillations [360] . 

In many neural population models, such as the well known Wilson-Cowan model [349] . it is 
assumed that the interactions are mediated by firing rates rather than action potentials (spikes) 
per se. To see how this might arise we rewrite Q in the equivalent form using a sum of Dirac 
d-functions 



Identifying the right hand side of Q as a train of pre-synaptic spikes motivates the form of a 
phenomenological rate model in the form 


Qs = /, (5) 

with / identified as a hring rate and Q identified as the differential operator (1 -|- Q!“^d/dt)(l -|- 
/3“^d/df). At the network level it is then common practice to close this system of equations by 
specifying / to be a function of pre-synaptic activity. A classic example is the Jansen-Rit model 
M, which describes a network of interacting pyramidal neurons (P), inhibitory interneurons (I) 
and excitatory interneurons (E), and has been used to model both normal and epileptic patterns 
of cortical activity [326] . This can be written in the form 

QeSP = f{sE — Sl), QeSE = C2f{CiSp) + A, QlSj = Cif{CsSp), 

which is a realisation of the structure suggested by Q , with the choice 


1 + Q-r(v-vo) ’ 
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and Qa = {l + /3~^d/dt)'^f3a/Aa for a G {E,I}. Here A is an external input. When this is a constant 
we obtain the bifurcation diagram shown in Fig. Oscillations emerge via Hopf bifurcations and it 
is possible for a pair of stable periodic orbits to coexist. One of these has a frequency in the alpha 
band and the other is characterised by a lower frequency and higher amplitude. Recently a network 
of such modules, operating in the alpha range and with additive noise, has been used to investigate 
mechanisms of cross-frequency coupling between brain areas USD]. Neural mass models have also 
previously been used to model brain resonance phenomena [300] , for modelling of epileptic seizures 
|45l 13191 165j . and are very popular in the Neuroimaging community for model driven EEG/fMRI 
fusion [330]. 



Eigure 6: Bifurcation diagram for the Jansen-Rit model with /3e = 100 s“^, /3/ = 50 s“^, Ae = 3.25 
mV, Aj = 22 mV, = 5 s“^, vq = 6 mV, r = 0.56 mV“^, Ci = 135, C 2 = O. 8 C 1 , C 3 = 0.25Ci = C 4 . 
Solid red (black) lines represent stable (unstable) fixed points. Green (blue) points denote the 
amplitude of stable (unstable) periodic orbits that emerge via Hopf bifurcations. The inset shows 
the co-existence of two stable periodic orbits at A = 125 Hz. 


Now that we have introduced some oscillator models for neurons and neural populations it is 
appropriate to consider the set of tools for analysing their behaviour at the network level. 

3 Dynamical systems approaches to collective behaviour 

We give a brief overview of some dynamical systems approaches, concepts and techniques that 
can be used to understand collective behaviour that spontaneously appears in coupled dynamical 
system models used for neuroscience modelling. We do not give a complete review of this area but 
try to highlight some of the approaches and how they interact; some examples of applications of 
these to neural systems are given in later chapters. 

3.1 Synchrony and asynchrony 

A set of N coupled nonlinear systems represented by ODEs can be written 

^Xi = gi{xi;xi,... ,xn), (6) 


13 













where Xi G represent the state space of the individual systems whose evolution gi is affected 
both by the current state of the system and by the states of those coupled to that system. In 
the mathematics literature, the Xi are often called “cells” though we note that the Xi may include 
degrees of freedom in the coupling as well as variables such as membrane potential that reflect the 
state of the “biological cell”. Note there is potential for notational confusion here: to clarify this 
we write 

X G Xi G G M. 

One of the most important observations concerning the collective dynamics of coupled nonlinear 
systems relates to whether the collection behaves as one or not - whether there is an attracting 
synchronous state, or whether more complex spatio-temporal patterns such as clustering appear. 
There is a very large literature, even restricting to the case of applications of synchrony, and one 
where we cannot hope to do the whole area justice. We refer in particular to |13L I255j . Various 
applications of synchrony of neural models are discussed, for example, in [52l [STJ [63l ESI 123711248L 
12561 12581 I276j while there is a large literature (e.g. [329]) discussing the role of synchrony in 

neural function. Other work looks for example at synchronisation of groups of networks [299] and 
indeed synchrony can be measured experimentally [270] in groups of neurons using dynamic patch 
clamping. 

We discuss some of the types of behaviour that can emerge in the collective dynamics and 
the response of partial synchronised states to external forcing. Clearly, any system of N coupled 
dynamical systems can be written in the form 

^Xj = fi{xi) -b egi{xi; xi,..., xn), (7) 


where each system is parametrised by x* G e = 0 corresponds to decoupling the systems and 
the functions gi represent drive from other systems on the zth system. Many approaches start with 
the less general case 


d 



N 

fi{xi) + e'^WijGij{xi,Xj), 
j=0 


( 8 ) 


which can be justified for system where the collective interaction between systems can be broken 
into “pairwise” interactions Gij that are summed according to some linear weights Wij (some of 
which may be zero) that represent the strength of the couplings and e the strength of coupling 
of the network. Note that there is clearly no unique way to write the system in this form; more 
specifically, one can without loss of generality choose e = 1, Wij = 1 by suitable choice of Gij. On 
the other hand it can be useful to be able to e.g. modulate the strength of the coupling across the 
whole network independently of changing individual coupling strengths. Similarly it is often useful 
to specialise to a case such as Gij = G and have similar interactions. 


3.2 Networks, motifs and coupled cell dynamics 

We focus now on the dynamics of pairwise coupled networks such as Q as this form is assumed 
in most cases. Under an additional assumption that the coupling between the oscillators is of the 
same type and either present or absent, one can consider uniformly weighted coupling of the form 

^ij — 9-^ij ) 
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where g is fixed and Aij is an adjacency matrix of the graph of interactions, i.e. Aij = 1 if there is 
a link from j to i and 0 otherwise, or more generally 


^ij — 9ij -^ij ) 


where gtj > 0 represents the strength of coupling and Aij the adjacency matrix. It is clear that the 
coupling structure as represented in Aij will influence the possible dynamics on the network and 
to make further progress it is useful to restrict to particular network structures. Some important 
literature on network structures is reviewed for example in |234j . while |13j reviews work on syn¬ 
chrony in complex networks up to the time of their publication; recent work includes for example 
|92j . For a recent review of the application of complex networks to neurological disorders in the 
brain, see [303]. 

It is interesting to try and understand the effect of network structure on synchrony, so we 
briefly outline some basic measures of network structure. The in-degree of the node i is the number 
of incoming connections (i.e. = Ylj^ji)^ while the out-degree is the number of outgoing 

connections (i.e. ciout(*) = Ylj ^ji) distribution of these degrees is often used to characterise 

a large graph. A scale-free network is a large network where the distribution of in (or out) degrees 
scales as a polynomial. This can be contrasted with highly structured homogeneous networks (for 
example on a lattice) where the degree may be the same at each node. Other properties commonly 
examined include the clustering properties and path lengths within the graph. There are also 
various measures of centrality that help one to determine the most important nodes in a graph - 
for example the betweenness centrality is a measure of centrality that is the probability that a given 
node is on the shortest path between two uniformly randomly chosen nodes [234] . As expected, the 
more central nodes are typically most important if one wishes to achieve synchrony in a network. 

Other basic topological properties of networks that are relevant to their dynamics include, for 
example, the following, most of which are mentioned in [131 1234] : The network is undirected if 
Aij = Aji for all i,j, otherwise it is directed. We say nodes j and i in the network Aij are path- 
connected if for some n there is a path from j to i, i.e. {A^)ij 0 for some n. The network is 
strongly connected if for each i,j is path connected in both directions while it is weakly connected 
if we replace Aij by max(Ajj, Aji) (i.e. we make the network undirected) and the latter network is 
strongly connected. N.B. There may be some confusion in that strong and weak connectivity are 
properties of a directed network - while strong and weak coupling are properties of the coupling 
strengths for a given network! The diameter of a network is the maximal length of a shortest 
path between two points on varying the endpoints. Other properties of the adjacency matrix 
are discussed for example in |34j where spectral properties of graph Laplacians are linked to the 
problem of determining stability of synchronised states. Other work we mention is that of Pecora 
et al. [25211250] on synchronisation in coupled oscillator arrays (and see § 4.1), while |274j explores 


4.2). 


the recurrent appearance of synchrony in networks of pulse-coupled oscillators (and see 

Finally, we mention network motifs - these are subgraphs that are “more prevalent” than others 
within some class of graphs. More precisely, given a network one can look at the frequency with 
which a small subgraph appears relative to some standard class of graphs (for example random 
graphs) and if a certain subgraph appears more often than expected, this characterises an important 
property of the graph |219| . Such analysis has been used in systems biology (such as transcription 
or protein interaction networks) and has been applied to study the structure in neural systems (see 
for example |3n2[ 1301] ) and the implications of this for the dynamics. The have also been used to 
organise the analysis of the dynamics of small assemblies of coupled cells, see for example [1511139] . 
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3.3 Weak and strong conpling 

Continuing with systems of the form Q or Q, if the coupling parameter e is, in some sense, small 
we refer to the system as “weakly coupled”. Mathematically, the weak coupling approximation is 
very helpful because it allows one to use various types of perturbation theory to investigate the 
dynamics |142j . For coupling of limit cycle oscillators it allows one to greatly reduce the dimension 
of phase space. Nonetheless, many dynamical effects (e.g. “oscillator death” where the oscillations 
in one or more oscillators are completely suppressed by the action of the network |103| ) cannot occur 
in the weak coupling limit and moreover real biological systems that often have “strong coupling”. 
We will return to this topic to discuss oscillator behaviour. One can sometimes use additional 
structure such as weak dissipation and weak coupling of the oscillators to perform a semi-analytic 
reduction to phase oscillators; see for example [THl I283j . 


3.4 Clusters, exact and generalised synchrony 

If one has a notion of synchrony between the systems of Q, it is possible to discuss clustering 
according to mutual synchrony. Caution needs to be exercised whenever discussing synchrony - 
there are many distinct notions of synchrony which may be appropriate in different contexts and, 
in particular, synchrony is typically a property of a particular solution at a particular point in time 
rather than a property of the system as a whole. 

An important case of synchrony is exact synchrony, we say Xi{t) and Xj{t) are exactly synchro¬ 
nised if Xi{t) = Xj{t) for all t. Generalised synchrony is, as the name suggests, much more general 
and corresponds to there simply being a functional relationship of the form Xi{t) = F{xj{t)). For 
oscillators, phases can be used to define additional notions such as phase and frequency synchrony: 


see 


6.1 


Although we focus mostly on individual units with simple (periodic) dynamics, if the units have 
more complex dynamics (such as “chaotic oscillators”) one can understand synchrony of the cells 
by analysis of the linearised differences between oscillators, and there is a sizeable literature on 
this; see the review |43) . or [23 Ij for clusters in a system of globally coupled bistable oscillators. In 
the case of two linearly coupled identical chaotic oscillators 

= /(xi) -h e(x2 - xi), ^X2 = /(x2) + e(xi - X2), 


where (xi, X2) G if the individual oscillator dx/dt = /(x) has a chaotic attractor A C then 
the coupled system will have an attracting exactly synchronised attractor A = {(x,x) : x G A} 
only if the coupling e is sufficiently large in relation to the maximal Liapunov exponent of the 
synchronous state [43] . 


3.5 Synchrony, dynamics and time delay 


An area of intense interest is the role of time delay in the collective dynamics of coupled systems. In 
the neural context it is natural to include propagation delay between neurons explicity, for example 
in models such as 


N 


^Xi(t) = fi{xi{t)) + e'^WijGij{xi{t),Xj{t - r)), 
i=i 
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where the delay time r > 0 represents the time of propagation of the signal from one neuron to 
another. This presents a number of serious mathematical challenges, not least due to the fact that 
delay dynamical systems are infinite dimensional: one must specify the initial condition over a time 
interval t € [to — r, to] in order to have any chance of uniquely defining the future dynamics; this 
means that for a state to be linearly stable, an infinite number of eigenvalues need to have real part 
less than zero. 

Nonetheless, much can be learned about stability, control and bifurcation of dynamically syn¬ 
chronous states in the presence of delay; for example [HI 11381 IHTl 125811191j . and the volume [35] 
includes a number of contributions by authors working in this area. There are also well-developed 
numerical tools such as DDE-BIFTOOL |94l 1291] that allow continuation, stability and bifurcation 
analysis of coupled systems with delays. For an application of these techniques to the study of a 
Wilson-Cowan neural population model with two delays we refer the reader to m- 

3.6 A short introduction to symmetric dynamics 

Although no system is ever truly symmetric, in practise many models have a high degree of sym- 
metryQ Indeed many real world networks that have grown (e.g. giving rise to tree-like structures) 
are expected to be well approximated by models that have large symmetry groups [205]. 

Symmetric (more precisely, equivariant) dynamics provides a number of powerful mathematical 
tools that one can use to understand emergent properties of systems of the form 

= f{x), (9) 

with X G . We say (§ is equivariant under the action of a group T if and only if f{gx) = gf{x) 
for any 5 G T and x G . There is a well developed theory of dynamics with symmetry; in 
particular see [1211112511122j . These give methods that help in a number of ways: 

• Description: one can identify symmetries of networks and dynamic states to help classify 
and differentiate between them 

• Bifurcation: there is a well-developed theory of bifurcation with symmetry to help un¬ 
derstand the emergence of dynamically interesting (symmetry broken) states from higher 
symmetry states 

• Stability: bifurcation with symmetry often gives predictions about possible bifurcation sce¬ 
narios that includes information about stability 

• Generic dynamics: symmetries and invariant subspaces can provide a powerful structure 
with which one can understand more complex attractors such as heteroclinic cycles 

• Design: one can use symmetries to systematically build models and test hypotheses 

The types of symmetries that are often most relevant for mathematical modelling of finite networks 
of neurons are the permutation groups, i.e. the symmetric groups and their subgroups. Nonetheless, 
continuum models of neural systems may have continuous symmetries that influence the dynamics 
and can be used as a tool to understand the dynamics; see for example m- 

^Indeed, the human brain consists of the order of 10^^ neurons, but of the order of 100 — 1000 types 
http://neuromorpho.org/neuroMorpho/index.jsp meaning there is a very high replication of cells that are only 
different by their location and exact morphology. 
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Name 

Symbol 

Comments 

Full permutation 

Sn 

Global or all-to-all coupling [23l [5l] 

Undirected ring 

Bn 

Dihedral symmetry [231 IM] 

Directed ring 

Zn 

Cyclic symmetry [23l [5l| 

Polyhedral networks 

G 


Lattice networks 

Gi X Ga 

Gi and G 2 could be Dfc or 

Hierarchical networks 

G 1 IG 2 

Gi is the local symmetry, G 2 the 
global symmetry, and 1 is the wreath 
product [9T|. 


Table 1: Some permutation symmetry groups that have been considered as examples of symmetries 
of coupled oscillator networks. 

3.7 Permutation symmetries and oscillator networks 

We review some aspects of the equivariant dynamics that have proven useful in coupled systems that 
are relevant to neural dynamics - see for example [231 08] . In doing so we mostly discuss dynamics 
that respects some symmetry group of permutations of the systems. The full permutation symmetry 
group (or simply, the symmetric group) on N objects, Sn, is defined to be the set of all possible 
permutations of N objects. Formally it is the set of permutations a : {1,..., A^} —)• {1,... ,N} 
(invertible maps of this set). To determine effects of the symmetry, not only the group must be 
known but also its action on phase space. If this action is linear then it is a representation of 
the group. The representation of the symmetry group is critical to the structure of the stability, 
bifurcations and generic dynamics that are equivariant with the symmetry. 

For example, if each system is characterised by a single real variable, one can view the action 
of the permutations on as a permutation matrix 

o-Jij I g otherwise ’ 

for each ct G F; note that M„Mp = M^p for any cr, p G F. Table lists some commonly considered 
examples of symmetry groups used in coupled oscillator network models. 

More generally, for ([^ equivariance under an action of F means that for all cr G F, x G 
and i = 1,..., iV we have 

fa(i) (^a-(i )) T ) • • • ) ®w) fi (®cr(i)) T ^9i (®cr(i)) ^a{l )) • ■ • ) ^a{N )) • 

A simple consequence of this is: if F acts transitively on {1,..., (i.e. if for any i and j there is a 

cr G F such that a{j) = i) then all oscillators are identical, i.e. fi{xi) = F{xi) for some function F. 

The presence of symmetries means that solutions can be grouped together into families - given 
any x the set Fx := G F} is the group orbit of x and all points on this group orbit will 

behave in dynamically the same way. 

3.8 Invariant subspaces and symmetries 

Given a point x G we define the isotropy subgroup (or simply the symmetry) of x under the 
action of F on to be 

Sx := e T : gx = x}. 
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Figure 7: Suppose the six-element group F = D 3 of symmetries of the equilateral triangle acts on 
generated by a rotation g 2 and a reflection gi in the x-axis. The group orbit of the point u 
that is not fixed by any symmetries also has six elements (shown by filled circles), while any group 
orbit of a point v that is fixed by a symmetry (e.g. gi) has correspondingly fewer points (shown 
by open circles) in the group orbit. Bifurcation of equilibria with more symmetry typically leads 
to several equilibria with less (“broken”) symmetry. 


This is a subgroup of F, and the set of these groups forms a lattice (the isotropy lattice) by inclusion 
of subgroups. They are dynamically important in that for any trajectory x{t) we have Sa,(o) = ^x(t) 
for all t. A converse problem is to characterise the set of all points with a certain symmetry. If H 
is a subgroup (or more generally, a subset) of F then the fixed point space of H is defined to be 

Fix(77) := {x G : gx = x for all (7 G 77}. 

Because all x G Fix(77) have symmetry H these subspaces are dynamically invariant. See Fig. 
that illustrates this principle. Typical points x G Fix(77) have = 77, however some points may 
have more symmetry; more precisely, if 77 C 77 are isotropy subgroups then Fix(77) D Fix(77); and 
the partial ordering of the isotropy subgroups corresponds to a partial ordering of the fixed point 
subspaces. 

One can identify similar types of isotropy subgroup as those that are conjugate in the group, 
i.e. we say 77 i and H 2 are conjugate subgroups of F if there is a 5 G F such that gHi = H 2 g. If 
this is the case, note that 


gFix(77i) 


g{x G : hx = x for all h G 77 i} 

{x G : hg~^x = g~^x for all h G 77 i} 

{x G : ghg~^x = x for all h £ Hi} 

{x G : hx = X for all h G 772} = Fix(772), 


meaning that the fixed point spaces of conjugate subgroups (and the dynamics on them) are in 
some sense equal. 

By identifying symmetries up to conjugacy allows for a considerable reduction of the number 
of cases one nes to consider; note that conjugate subgroups must have fixed point subspaces of the 
same dimension where essentially the same dynamics will occur. 

The fixed point subspaces are often used (implicitly or explicitly) to enable one to reduce the 
dimension of the system and to make it more analysable. As an example, to determine the existence 
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of a fully exact synchrony solution one only need to suppose Xi{t) = x{t) and determine whether 
there is such a solution x{t) for the system Q. 

For periodic orbits, the symmetry of points on the orbit to symmetries of the orbit as an 
invariant set are as follows. Suppose P is a periodic orbit (which can be viewed as a “loop” in 
phase space Let K denote the symmetries that fix all points on P (the “spatial symmetries”) 
and H denote the symmetries that fix P as a set (the “spatio-temporal symmetries”); note that K 
will be a subgroup of H. Finally, let 


Lk ■= y Fix({5f}) nFix(P:), 

g^K 

be the set of points in phase space that have strictly more symmetry than K. 

Theorem 3.1 \12^ Theorem 3-4] Consider ODEs on M"' with a given finite symmetry group F. 
There is a periodic orbit P with spatial symmetries K and spatio-temporal symmetries H if and 
only if 

• H/K is cyclic 

• K is an isotropy subgroup 

• dimFix(iF) > 2 

• H fixes a connected component o/Fix(iF)/Lii', where Lk is defined as above. 


One way of saying this is that the only possible spatio-temporal symmetries of periodic orbits 
are cyclic extensions of isotropy subgroups. Further theory, outlined in |122j . shows that one can 
characterise possible symmetries of chaotic attractors; these may include a much wider range of 
spatio-temporal symmetries (P, K) including some that do not satisfy the hypotheses of Theo¬ 
rem This means that the symmetries of attractors may contain dynamical information about 
the attractor. 

3.9 Bifurcations with symmetry and genericity 

Bifurcations of ODEs can be classified and analysed by codimension according to methods for 
example in texts |180(I12f^ . This involves the following steps: 

(a) Identification of the marginally unstable modes (the directions that are losing stability). 

(b) Reduction to a centre manifold parametrised by the marginally unstable modes (generically 
this is one or two dimensional when only one parameter is varied). 

(c) Study of the dynamics of the normal form for the bifurcation under generic assumptions on 
the normal form coefficients. 

Indeed, the only generic codimension one local bifurcations (i.e. the only one-parameter bifurcations 
of equilibria that will not split into a number of simpler bifurcations for some small perturbation on 
the system) are the saddle-node (also called fold, or limit point) and the Hopf (also called Andronov- 
Hopf) bifurcation. Additional more complicated bifurcations can appear at higher codimension. 
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This classification by codimension has enabled development of a powerful set of numerical tools 
to help analysis of such systems, not just for local bifurcations of equilibria but also some global 
bifurcations (in particular, periodic orbit and homoclinic bifurcations). Particular packages to do 
this include AUTO [169] . MatCont [8^, CONTENT [179j and XPPAUT [99] (which includes an 
implementation of AUTO). 

If we restrict to symmetry preserving perturbations, a much wider range of bifurcations can 
appear at low codimension - this is because the symmetry can cause a marginal mode of instability 
in one direction to appear simultaneously in many other directions meaning that 

(a’) Identification of the marginally unstable modes (symmetry means there can generally be 
several of these that will go unstable at the same time). 

(b’) Reduction to a centre manifold parametrised by the marginally unstable modes (these are 
preserved by the action of the symmetries and may be of dimension greater than two even 
for one parameter bifurcations). 

(c’) Study of the dynamics of the normal form for the symmetric bifurcation under generic as¬ 
sumptions on the normal form coefficients (the symmetries mean that some coefficients may 
be zero, some are constrained to be equal while others may be forced to satisfy nontrivial and 
sometimes obscure algebraic relationships). 

These factors conspire to make symmetric bifurcations rich and interesting in behaviour - even 
in codimension one it is possible for heteroclinic cycles or chaos to bifurcate from high symmetry 
solutions. However, the same factors mean that many features cannot be caught by numerical 
path-following packages such as those listed above - the degeneracies mean that many branches 
may coming out of the bifurcation; it is generally a challenge to identify all of these. Essentially, 
bifurcation theory needs to be developed in the context of the particular group action. Examples 
of some consequences of this for weakly coupled oscillator networks with symmetries are considered 
in § [6] 

3.10 Robust heteroclinic attractors, cycles and networks 

The presence of symmetries in a dynamical system can cause highly nontrivial dynamics even away 
from bifurcation points. Of particular interest are robust invariant sets that consist of networks of 
equilibria (or periodic orbits, or more general invariant sets) connected via heteroclinic connections 
that are preserved under small enough perturbations that respect the symmetries [ITOj . These 
structures may be cycles or more generally networks. They can be robust to perturbations that 
preserve the symmetries and indeed they can be attracting [TTniT^ . We are particularly interested 
in the attracting case in which case we call these invariant sets heteroclinic attractors and trajectories 
approaching such attractors show a typical intermittent behaviour - periods that are close to the 
dynamics of an unstable saddle-type invariant set, and switches between different behaviours. 

In higher dimensional systems, heteroclinic attractors may have subtle structures such as “depth 
two connections” |28|, “cycling chaos” where there are connections between chaotic saddles |85| 
11221 [31] and “winnerless competition” |263L 1261] . Related dynamical structures are found in the 
literature in attractors that show “chaotic itinerancy” or “slow switching”. Such complex attractors 
can readily appear in neural oscillator models in the presence of symmetries and have been used to 
model various dynamics that contribute to the function of neural systems; we consider this, along 
with some examples, in § [7| 
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3.11 Groupoid and related formalisms 

Some less restrictive structures found in some coupled dynamical networks also have many of 
the features of symmetric networks (including invariant subspaces, bifurcations that appear to be 
degenerate, and heteroclinic attractors) but without necessarily having the symmetries. 

One approach |307] has been to use a structure of groupoids - these are mathematical structures 
that satisfy some, but not all, of the axioms of a group and can be useful in understanding the 
constraints on the dynamics of coupled cell systems of the form Q. A groupoid is similar to a 
group except that the composition of two elements in a groupoid is not always defined, and the 
inverse of a groupoid element may only be locally defined. This formalism can be used to describe 
the permutations of inputs of cells as in |3n7( 1123] . 

Suppose that we have Q with cells C = {1,..., A"} and suppose that there are connections T, 
i.e. are pairs of cells (f,j) in C such that cell i appears in the argument of the dynamics of cell j. 
We say 

/(j) = {f G C : (i, j) G T}, 

is the input set of cell j and there is a natural equivalence relation ~7 defined by j k if there is 
a bijection (an input automorphism) 


with j3{j) = k such that for all i G I{j) we have {i,j) ~_e {I3{i),k). The set B{j,k) of input 
automorphisms of this type and the set of all such input automorphisms has the structure of a 
groupoid [Mi- 

Given a coupling structure of this type, an admissible veetor field is a vector field on the 
product space of all cells that respects the coupling structure, and this generalises the idea of an 
equivariant vector held in the presence of a symmetry group acting on the set of cells. The dynamical 
consequences of this have a similar havour to the consequences one can hnd in symmetric systems 
except that fewer cases have been worked out in detail, and there are many open questions. 

To illustrate, consider the system of three cells 

5(xi,X2,X3), 

g{x2,xi,xfi), 
h{x3,xi), 

where g{x,y,z) = g{x,z,y)] this is discussed in [1231 §5] in some detail. This has a coupling 
structure as shown in Fig. In spite of there being no exact symmetries in the system there is a 
nontrivial invariant subspace where xi = X 2 - In the approach of m the dynamically invariant 
subspaces that can be understood in terms of the balanced colourings of the graph where the cells 
are grouped in such a way that the inputs are respected - this corresponds to an admissible pattern 
of synchrony. 

The invariant subspaces that are forced to exist by this form of coupling structure have been 
called polydiagonals in this formalism, which correspond to clustering of the states. For every 
polydiagonal one can associate a quotient network by identifying cells that are synchronised, to 
give a smaller network. As in the symmetric case the existence of an invariant subspace does 


dt"^ = 


dt"^ = 


dt 


X3 = 


22 











Figure 8: Left: a system of three coupled cells with two cell types (indicated by the circle and 
square) coupled in a way that there is no permutation symmetry but there is an invariant subspace 
corresponding to cells 1 and 2 being synchronised. The different line styles show coupling types 
that can potentially be permuted (after Fig. 3]). Middle: the quotient two-cell network 

corresponding to cells 1 and 2 being synchronised. Right: The same network structure shown using 
the notation of [9] . 


not guarantee that it contains any attracting solutions. Some work has been done to understand 
generic symmetry breaking bifurcations in such networks - see for example m, or spatially periodic 
patterns in lattice networks [HS] . Variants of this formalism have been developed to enable different 
coupling types between the same cells to be included. 

Periodic orbits in such networks can also have interesting structures associated with the presence 
of invariant subspaces. The so-called rigid phase conjecture [30811123] . recently proved in [120) . 
states that if there is a periodic orbit in the network such that two cells have a rigid phase relation 
between (i.e. one that is preserved for all small enough structure-preserving perturbations) then 
this must be forced by either a symmetric perturbation of the cells in the network, or in some 
quotient network. 

An alternative formalism for discussing asymmetric coupled cell networks has been developed 
in [ml 0 0 m that also allows one to identify invariant subspaces. Each cell has one output and 
several inputs that may be of different types. These papers concentrate on the questions: (a) When 
are two cell networks formally equivalent (i.e. when can the dynamics of one cell network be found 
in the other, under suitable choice of cell)? (b) How can one construct larger coupled cell systems 
with desired properties by “inflating” a smaller system S, such that the larger system has S' as a 
quotient? (c) What robust heteroclinic attractors exist in such systems? 


4 Coupled limit cycle oscillators 


In addition to variants on the systems in § 0 we mention that some nonlinear “conceptual” planar 
limit cycle oscillators are also studied in neural models. These include: 




/ 2 \ ^ 2 
X + (ao -|- aix ) —X -|- uj 




X -t- (ao -|- aix‘^)—x -|- -|- a2X^ = 0 

—z = (A -|- iu})z -|- (ao -|- iai)\z\^z 


van der Pol, 

van der Pol-Dufhng, 

Stuart-Landau oscillator, 


where x is real, 2 is complex and the coefficients A, u and aj are real constants. 

If the coupling between two or more limit cycle oscillators is relatively large, it can affect not 
only the phases but also the amplitudes, and a general theory of strongly interacting oscillators 
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is likely to be no more or less complicated than a general theory of general nonlinear systems. 
However, the theory of weak-coupling is relatively well developed (see § and § and specific 
progress can sometimes be made for the special choice of neuron model. Examples where one can 
do this include IF (see § 4.3), piece-wise linear models such as McKean caricatures of FHN, ML 


m and singularly perturbed relaxation oscillators with linear m or fast threshold modulation 
coupling |298j . 

For linear coupling of planar oscillators, much is known about the general case |15l I177| . If 
the linear coupling is proportional to the difference between two state variables this is referred to 
as “diffusive”, and otherwise it is called “direct”. The difference between the two cases is most 
strongly manifest when considering the mechanism of oscillator death (see § 3.3). The diffusive 


case is more natural in a neuroscience context as it can be used to model electrical gap junction 
coupling (which depends on voltage-differences). The existence of synchronous states in networks 
of identical units is inherited from the properties of the underlying single neuron model since in 
this case coupling vanishes, though the stability of this solution will depend upon the pattern of 
gap-junction connectivity. 

Gap junctions are primarily believed to promote synchrony, though this is not always the case 
and they can also lead to robust stable asynchronous states IM], as well as “bursting” generated by 
cyclic transitions between coherent and incoherent network states US]. For work on gap junctions 
and their role in determining network dynamics see for example |260L I254L I%1 [8^ 12081 172] . 


4.1 Stability of the synchronised state for complex networks of identical systems 


There is one technique specific to the analysis of the synchronous state in a quite large class of 
network models that is valid for strongly coupled identical systems, namely the master stability 
function (MSF) approach. For networks of coupled systems or oscillators with identical compo¬ 
nents the MSF approach of Pecora and Carroll [253] can be used to determine the stability of 
the synchronous state in terms of the eigenstructure of the network connectivity matrix. To in¬ 
troduce the MSF formalism it is convenient to consider N nodes (oscillators)!^ and let Xi G 
be the d dimensional vector of dynamical variables of the ith node with isolated (uncoupled) dy¬ 
namics dxi/dt = F{xi), with i = The output for each node is described by a vector 

function 77 G (which is not necessarily linear). For example for a three dimensional system with 
X = {x^^\ x^'^\ x^^'i) and linear coupling only occurring through the -component then we would 
set H{x) = (0,0,x*^^)). For a given adjacency matrix Aij and associated set of connection strengths 
9ij > 0 and a global coupling strength a the network dynamics of coupled identical systems, to 
which the MSF formalism applies, is specified by 


d ^ 

—Xi = F{xi)+aY^ 

i=i 


N 


AijQij [H{xj) - H{xi)] = F{xi) - a'^QijH{xj). 


Here the matrix Q with blocks Gij has the graph-Laplacian structure Gij = —Aijgij + 6ij ^ikdik- 
The N — 1 constraints xi{t) = X 2 {t) = ... = X]sf{t) = s{t) define the (invariant) synchronisation 
manifold, with s{t) a solution in of the uncoupled system, namely ds/dt = F{s). Although we 
will not discuss in detail here, we assume that the asymptotic behaviour s{t) is such that averages 
along trajectories converge, i.e. the behaviour of s{t) for the uncoupled system is governed by a 
natural ergodic (SRB) measure for the dynamics. 

^In this section we assume little about the dynamics of the nodes - they may be “chaotic oscillators”. 
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To assess the stability of this state we perform a linear stability analysis expanding a solution 
as Xi{t) = s{t) + 5xi{t) to obtain the variational equation 


N 


(5—Xj = DF{s)5xi — aDH{s) Gijdxj. 

Clt 

J=1 


Here DF(s) and DH{s) denote the Jacobian of F{s) and H{s) around the synchronous solution 
respectively. The variational equation has a block form that can be simplified by projecting 5x 
into the eigenspace spanned by the (right) eigenvectors of the matrix Q. This yields a set of N 
decoupled equations in the block form 

= [DF{s) - a\iDH{s)\ iu l = 

where is the Ith (right) eigenmode associated with the eigenvalue A/ of G (and DF{s) and 
DH{s) are independent of the block label 1). Since — 0 there is always a zero eigenvalue, 

say Ai = 0, with corresponding eigenvector (1,1,, 1), describing a perturbation parallel to the 
synchronisation manifold. The other N — 1 transverse eigenmodes must damp out for synchrony 
to be stable. For a general matrix Q the eigenvalues Xi may be complex, which brings us to 
consideration of the system 

^^=[DF{s)-aDH{s)]^, a = aXiGC. (10) 


For given s{t), the MSF is defined as the function which maps the complex number a to the greatest 
Lyapunov exponent of (10). The synchronous state of the system of coupled oscillators is stable if 
the MSF is negative at a = aXi where A; ranges over the eigenvalues of the matrix Q (excluding 
Ai =0). 

For a ring of identical (or near identical) coupled periodic oscillators in which the connections 
have randomly heterogeneous strength Restrepo et al. |269j have used the MSF method to deter¬ 
mine the possible patterns at the desynchronisation transition that occurs as the coupling strengths 
are increased. Interestingly they demonstrate Anderson localisation of the modes of instability, and 
show that this could organise waves of desynchronisation that would spread to the whole network. 
For a further discussion about the use of the MSF formalism in the analysis of synchronisation of 
oscillators on complex networks we refer the reader to |13( 1259) , and for the use of this formalism 
in a non-smooth setting see [320) . This approach has recently been extended to cover the case 
of cluster states by making extensive use of tools from computational group theory to determine 
admissible patterns of synchrony [251] (and see also 


3.11) in unweighted networks. 


4.2 Pulse-coupled oscillators 

Another example of a situation in which analysis of network dynamics can be carried out without 
the need for any reduction or assumption is that of pulse coupled oscillators, in which interactions 
between neurons are mediated by instantaneous “kicks” of the voltage variable. 

Networks of N identical oscillators with global (all-to-all) strong pulse coupling were first studied 
by Mirollo and Strogatz |222] . They assumed that each oscillator is dehned by a state variable v 
and is of integrate-and-fire type with threshold vth = 1 and reset value u/j = 0. When oscillator 
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i in the network fires the instantaneous pulsatile coupling pulls all other oscillators j 7 ^: i up by a 
hxed amount e or to firing, whichever is less, i.e. 


If Vi{t) = 1 then = min(l, Vj{t) + e) for all j / i. 


Mirollo and Strogatz assume that the coupling is excitatory (e > 0). If m oscillators hre simulta¬ 
neously then the remaining N — m oscillators are pulled up by me, or to firing threshold. 

In the absence of coupling each oscillator has period A and there is a natural phase variable 
= t/A mod 1 such that (j) = 0 when u = 0 and cj) = 1 when v = 1. Mirollo and Strogatz 
further assume that the dynamics of each (uncoupled) oscillator is governed by v{t) = f{4>) where 
/ is a smooth function satisfying /(O) = 0, /(I) = 1, f'{4>) > 0 and /"(</>) < 0 for all cj) G [0,1]. 
Because of these hypotheses on /, it is invertible with inverse (p = g{v). 

Note that for the leaky (linear) integrate-and-fire model (LIF) where 


T—v = —V + I 
at 


we have f{(p) = /(I — where A = rln(//(/ — 1)), for / > 1, which satisfies the above 

conditions. However, quadratic IF models do not satisfy the concavity assumption. 

If an oscillator is pulled up to firing threshold due to the coupling and firing of a group of m 
oscillators which have already synchronised then the oscillator is ‘absorbed’ into the group and 
remains synchronised with the group for all time. (Here synchrony means firing at the same time). 
Since there are now more oscillators in the synchronised group, the effect of the coupling on the 
remaining oscillators is increased and this acts to rapidly pull more oscillators into synchronisation. 
Mirollo and Strogatz |222j proved that for pulsatile coupling and / satisfying the conditions above, 
the set of initial conditions for which the oscillators do not all become synchronised has zero measure. 
Here we briefly outline the proof for two pulse coupled oscillators. See Mirollo and Strogatz [222] 
for the generalisation of this proof to populations of size N. 

Consider two oscillators labelled A and B with 4>a and va denoting respectively the phase and 
state of oscillator A and similarly for oscillator B. Suppose that oscillator A has just fired so that 
pA = 0 and (pB = p as in Fig. (a). The return map R{p) is the phase of B immediately after A 
next fires. It can be shown that the return map has a unique, repelling fixed-point: 

Oscillator B reaches threshold when pA = 1 — p and an instant later B hres and the pulsatile 
coupling makes va = min(l,e -|- /(I — </>)). If = 1 then the oscillators have synchronised. 
Assuming that va = ^ + /(I — p) < I then pA = g{e -|- /(I — p)) and after one firing the system has 
moved from {pA^ Pb) = (0, p) to {h{p), 0 ) where h{p) = g{e + /(I — p)) is the firing map (see Fig. 
(c)). The return map is given by a further application of the firing map so that R{p) = h{h{p)). The 
assumption that e+f{l—p) < 1 is satisfied when e G [0,1) and p G (<5,1) where 5 = l—g{l—e). Thus 
the domain of h is {6, 1) and the domain of R is Since h is monotonically decreasing, 

6 < h~^[6) for e < 1 and the interval is nonempty. Thus on the whole of [0,1) the return map is 
defined as 

f 1 P> 

R{P) = { h{h{P)) P€{5,h-H6)), 

[0 p < 6. 

Since the points 0 and 1 are identified, if (/> < <5 or (/> > h~^{6) then the two oscillators will become 
synchronised. 


26 






(a) 


(b) 


(c) 




Figure 9: A system of two oscillators governed hy v = f{4>), and interacting by pulse coupling, (a) 
The state of the system immediately after oscillator A has fired, (b) The state of the system just 
before oscillator B reaches the firing threshold, (c) The state of the system just after B has fired. 
B has jumped back to zero, and the state of A is now min(l, e + /(I — (p)). 


It can be shown that almost all initial conditions eventually become synchronised since (i) R 
has a unique fixed point cp G {6,h~^{6)) and (ii) this fixed point is unstable (i.e. \R'{(p)\ > 1). To 
see that R has a unique fixed point, observe that fixed points (p are roots of F{(p) = cp — h{(p). Now 
F{5) = (5 — 1 < 0 and F{h~^{5)) = h~^{5) — <5 > 0 so F has a root in (<5, h~^{5)) and this root is 
unique since F'{(p) = 1 — h'{(p) > 2. 

Extensions to the framework of Mirollo and Strogatz include the introduction of a time delay in 
the transmission of pulses and the consideration of inhibitory coupling. It has been observed that 
delays have a dramatic effect on the dynamics in the case of excitatory coupling. Considering first 
a pair of oscillators, Ernst et al. [lOSj demonstrate analytically that inhibitory coupling with delays 
gives stable in-phase synchronisation while for excitatory coupling, synchronisation with phase lag 
occurs. As the number of globally coupled oscillators increases, so does the number of attractors 
which can exist for both excitatory and inhibitory coupling. 

In the presence of delays many different cluster state attractors can coexist. The dynamics 
settle down onto a periodic orbit with clusters reaching threshold and sending pulses alternately 
|108l 11091 1322j . Under the addition of weak noise when the coupling is inhibitory, the dynamics 
stay near this periodic orbit indicating that all cluster state attractors are stable [322]. However, 
the collective behaviour shows a marked difference when the coupling is excitatory. In this case, 
weak noise is sufficient to drive the system away from the periodic orbit and results in persistent 
switching between unstable (Milnor) attractors. 

These dynamics are somewhat akin to heteroclinic switching and the relationship between 
networks of unstable attractors and robust heteroclinic cycles has been addressed by a number 
of authors [33 Ea csni- In particular, Broer et al. usa highlight a situation in which there 
is a bifurcation from a network of unstable attractors to a heteroclinic cycle within a network 
of pulse coupled oscillators with delays and inhibitory coupling. They note that the model used 
in previous work [lOSl 11091 1322| is locally noninvertible since the original phase of an oscillator 
cannot be recovered once it has received an input which takes it over threshold causing the phase 
to be reset. Kirst and Timme employ a framework in which a reset function R{C) = c(p, 

c G [0,1] is introduced which ensures that the flow becomes locally time invertible when c > 0. 
They demonstrate that for c = 0 (where the locally noninvertible dynamics are recovered), the 
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system has a pair of periodic orbits and A 2 which are unstable attractors enclosed by the basin 
of each other. When c > Q, Ai and A 2 are non-attracting saddles with a heteroclinic connection 
from Ai to ^ 2 - Furthermore, there is a continuous bifurcation from the network of two unstable 
attractors when c = 0 to a heteroclinic two cycle when c > 0. 

For an interesting dynamical systems perspective on the differences between “kick” synchroni¬ 
sation (in pulsatile coupled systems) and “diffusive” synchronisation j3U6| and the lack of mathe¬ 
matical work on the former problem see [213j . For example, restrictions on the dynamics of sym¬ 
metrically coupled systems of oscillators when the coupling is time-continuous can be circumvented 
for pulsatile coupling leading to more complex network dynamics |158) . 

In the real world of synaptic interactions, however, pulsatile kicks are probably the exception 
rather than the rule, and the biology of neurotransmitter release and uptake is better modelled 
with a distributed delay process, giving rise to a post synaptic potential with a finite rise and fall 
time. For spike-time event driven synaptic models, described in § 2.5, analysis at the network level 


is hard for a general conductance based model (given the usual expectation the single neuron model 
will be high dimensional and nonlinear), though far more tractable for LIF networks, especially 
when the focus is on phase-locked states |337L 133311332j . Indeed in this instance many results can be 
obtained in the strongly coupled regime [49| . without recourse to any approximation or reduction. 


4.3 Synaptic coupling in networks of IF neurons 

The instantaneous reaction of one neuron to the firing of another, as in the pulse-coupled neurons 
above, does not account for the role of synapses in the transmission of currents. Bressloff and 
Coombes [l9] consider a network of N identical, LIF neurons that interact via synapses by trans¬ 
mitting spike trains to one another. Let fj(t), the state of neuron i at time t, evolve according to 

^Vi = -Vi +Ii + Xi{t), (11) 

where li is a constant external bias and Xi{t) is the total synaptic current into the cell. As before, 
we supplement this with the reset condition that whenever Vi = 1 neuron i fires and is reset to 
Vi = 0. The synaptic current Xi{t) is generated by the arrival of spikes from other neurons j and 
can be taken to have the form 

N 

W(t) = ej;n;,, ^ J(t-Tf), 

j=l rnGh 


where ewij represents the weight of the connection from the jth neuron to the ith neuron with e 
characterising the overall strength of synaptic interactions, T™ denotes the sequence of firing times 
of the jth neuron and J{t) determines the course of postsynaptic response to a single spike. A 
biologically motivated choice for J{t) is 

J{t) = r]{t)Q{t), r]{t) = 


where 0(t) = 1 if t > 0 and zero otherwise. Here rj is an alpha function (see also § 2.5) and the 
maximum synaptic response occurs at a nonzero delay t = a~^. Note that in the limit of large 
inverse rise time a, J{t) approximates a delta function (more like pulse-coupling). 

Bressloff and Coombes |49j show that the behaviour of the network of oscillators differs depend¬ 


ing on the strength of the coupling. This is another instance in which information is lost through 
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making weak-coupling assumptions. To see this one may integrate (11) between T” and 
exploiting the linearity of the equations and solving with variation of parameters, to obtain a map 
of the firing times. Since the drive Xi{t) depends upon all previous firing events of the other neuron 
this is an implicit map that relates all the network firing events to one another. It is convenient to 
introduce the set of inter-spike intervals (ISIs) = T” — TJ", so that we may write the firing 

map in the convenient form 


Ii{l - e 


-a; 


n+l,n 


N 


1 + ^ ^ Fij{A, 
j=l mGZ 


ri.+l,n 


A n,7n 


) = o, 


( 12 ) 


where > T™ for all j, m, and 


px 

Fij{x,y) = ewije~'^ eV(s-k?/)d 6 

Jo 


Phase-locked solutions may be found, for an arbitrary coupling strength e, using the ansatz Tj” = 
(m — 4>j)A for some self-consistent ISI A and constant phases (j)j. Substitution into (12) yields 


N 


where 


l = /i(l-e ^) + e'^WijK{(l)j - (l)i), 
i=i 


/■A 

7^(0) = e“^ / e^P{s + (/>A)ds, P{t) = J{t — mA), 

Jo 


(13) 


m£7j 


and P{t) = P{t -I- A). Choosing one of the phases, say ^i, as a reference then (13) provides a set 
of N equations for the unknown period A and the remaining N — 1 relative phases (jij — (pi. 

In order to investigate the linear stability of phase-locked solutions of equation (13), we consider 
perturbations of the firing times which we write in the form T™ = (m — (pj)A + Sp. Linearisation 
of the firing map (12) gives an explicit map for these perturbations that can be written as 

N 

(9x * ^ dy 


j=\ mGZ 


n+l,n 


where the partial derivatives of Fij are evaluated at the phase-locked state (A”’ , 

(A, (n — m)A -|- (pj — pi)A). For solutions of the form 5™ = this reduces to 


N 




(14) 


where - 1 + WijP{{Pj - Pi)A), Hij{X) = WijGHPj - Pi)A, A) - 5ij Y.k WikG{{Pk - 

0i)A, 1), 

and 


G(t, A) = V A-™e-^ / eV' 

mez -^0 


(s -I- t -I- mA)ds. 
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One solution to equation (14) is A = 1 with 6i = S for all i = 1,..., A^. This reflects the invariance 


of the dynamics with respect to uniform phase shifts in the firing times. Thus the condition for 


linear stability of a phase-locked state is that all remaining solutions A of equation (14) are within 
the unit disc. For A — 1 ~ 0{e), and e small, we may expand (14) as 


N 


(A - l){Ii - l)<5i = + 0(e2), 


i=i 

where = [WijK'{4>j — cpi) — {4’k — 4>i)]/^ and we have used the 

result that G{(p,l) = K'{(f>)/A. Suppose for simplicity that Ij = I > 1 for all i, so that A = 
ln[I/(/ — 1)] -|- 0(e). Then the weak-coupling stability condition is that ReAp < 0, where Ap, 
p = 1,..., N, are the eigenvalues of the N x N matrix with elements 

As an explicit example let us consider the synchronous state {(pi = cp for all i). From (13) 


we see that this is guaranteed to exist for the choice kFij = 7 and Ij = I[1 — eK{0)'j] for 
some constant / > 1 which sets the period as A = ln[//(/ — 1)]. Using the result that K'{(p) = 
— AK{(p) + AP{(pA)/1 the spectral problem for the synchronous state then takes the form 

N 

[(A - 1)(/ - 1 + el 7 A:'( 0 )/A) + e^K'{^)/A] 6i = eG(0, A) 

i=i 

We may diagonalise this equation in terms of the eigenvalues of the weight matrix, denoted by Vp 
with p = 1 ,..., A (and we note that 7 is the eigenvalue with eigenvector ( 1 , 1 ,..., 1 ) corresponding 
to a uniform phase-shift). Looking for non-trivial solutions then gives us the set of spectral equations 
£p{\) = 0 , where 


£p{X) = (A - 1)(/ - 1 + e/ 7 A'( 0 )/A) + e 7 A'( 0 )/A - ei/pG(0, A). 


(15) 


We may use (15) to determine bifurcation points defined by |A| = 1. For sufficiently small e, 
solutions to (15) will either be in a neighbourhood of the real solution A = 1 or in a neighbourhood 


of one of the poles of G(0, A). Since the latter lie in the unit disc, the stability of the synchronous 
solution (for weak coupling) is determined by eAr'(0)[Rer'p — 7 ] <0. For strong coupling the 
characteristic equation has been solved numerically in |49] to illustrate the possibility of Hopf 
bifurcations (A = e*^, 0 7 ^ 0, 0 7 ^ tt) with increasing |e|, leading to oscillator death in a globally 
coupled inhibitory network for a sufficiently slow synapse and bursting behaviour in asymmetrically 


coupled networks. Bifurcation diagrams illustrating these phenomenon are shown in Fig. 10 To see 


how these phenomena can occur from a more analytical perspective it is useful to consider the 
Fourier representation J{t) = (27r)“^ dw J(a;)e*‘^*, where J(w) = o?I{a P iw)^, so that G may 

be easily evaluated with A = e^ as 




1 J{ijJn- izj A){WnP zj A){(P - a ^), 




1 + iu)n zj^ 


ujn = 27rn/A, 


(16) 


from which it is easy to see a pole at z = —aA. This suggests writing 2 ; in the form 2 ; = —a(l-|-Kp)A 
and expanding the spectral equation in powers of a to find a solution. For small a we hnd from 


(16) that G(0,e^) = —a(l + «^p)(l — e ^)/{KpA). Balancing terms of order a then gives Kp = 


e^'p(l — e ^)/(A^(/ — 1)), where we use the result that G(0, e*^) = O(a^). Thus for small a. 
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Figure 10: Left: Region of stability for a synchronised pair of identical IF neurons with inhibitory 
coupling and collective period T = ln2. The solid curve |e| = ec(a) denotes the boundary of the 
stability region. Crossing the boundary from below signals excitation of the linear mode (1, —1) 
leading to a stable state in which one neuron becomes quiescent (oscillator death). For a > ao the 
synchronous state is stable for all e. The dashed curve denotes corresponds to the eigenvalue with 
ly = 1. Right: Plot of critical coupling €c as a function of a for various network sizes N. The critical 
inverse rise-time q;o(-^) is seen to be a decreasing function of N with ao{N) —)• 0 as iV —)• oo. 


namely slow synaptic currents, we have that K'(0) = 0, so that a weak-coupling analysis would 
predict neutral stability (consistent with the notion that a set of IF neurons with common constant 
input would frequency lock with an arbitrary set of phases). However, our strong coupling analysis 
predicts that the synchronous solution will only be stable if Re < 0 with = [—1 ± Kp]aA. 
Introducing the firing rate function / = 1/A then can be written succinctly as 


z„ = 


- 1 ± 




aA. 


Thus for an asymmetric network (with at least one complex conjugate pair of eigenvalues) it may 
occur that as |e| is increased a pair of eigenvalues determined by z^ may cross the imaginary axis 
to the right hand complex plane signalling a discrete Hopf bifurcation in the firing times. For a 
symmetric network with real eigenvalues an instability may occur as some Kp G M increases through 
1, signalling a breakdown of frequency locking. The above results (valid for slow synapses) can also 
be obtained using a firing rate approach, as described in [l9] . 

The results above, albeit valid for strong coupling, are only valid for LIF networks. To obtain 
more general results for networks of limit cycle oscillators it is useful to consider a reduction to 
phase models. 


5 Reduction of limit cycle oscillators to phase-amplitude and phase 
models 

Consider a system of the form 


—x = f{x) + eg{x,t), x€ 


(17) 
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such that for e = 0 the system possesses a periodic orbit 


r = {u{t) : t G M}, 

with minimal period T > 0 (such that u{t) = u{t + T) for all t G M but u{t) / u{s) for 0 < s < T). 
We will assume that T is attracting and hyperbolic, i.e. linearly stable so there one zero Floquet 
exponent and the others with negative real part. Then we say that 0 is a limit cycle oscillator. 
We will reduce this to a description that involves a phase that lives on a topological circle that 
can be thought of an interval [0, T) with the topology that comes from identifying the ends of the 
interval: e and T — e are close for e small. This circle is sometimes called the one-torus T. 

There are a number of conventions used in the literature to represent this phase. We discuss 
these conventions before looking at general reduction methods. One convention is to define a phase 
'd modulo T such that 

I-' 

where '&{u{t)) = t interpreted modulo T and such that i?(uo) = 0 for some chosen point uq G T. 
Another convention is to define a phase 9 modulo 2tt such that 


d 

s' 


= w, 


where co = 2tt/T is the angular frequency of the oscillator. Again we pick a point uq G T and 
require that 9 = 0 at u = uq. 

In the remainder of the article we will use and 9 as an indication of the convention we are 
using for phase, in the special case T = 2tt both conventions coincide and we use 9. Typically 
one of these descriptions is more convenient than the others but all can in principle be adapted to 
any phase reduction. Table expresses some features of the conventions which commonly used for 
phase variables. 

We now review some techniques of reduction which can be employed to study the dynamics of 
( |17[ ) when e 7 ^ 0 so that the perturbations may take the dynamics away from the limit cycle. In 
doing so we will reduce for example to an ODE for 'd{t) taken modulo T. Clearly any solution of 
an ODE must be continuous in t and typically tt^t) will be unbounded in t growing at a rate that 
corresponds to the frequency of the oscillator. Strictly speaking, the coordinate we are referring to 
in this case is on the lift of the circle T to a covering space M, and for any phase G [0, T) there 
are infinitely many lifts to M given by ?? -|- kT for k G Z. However, in common with most literature 
in this area we will not make a notational difference between whether the phase is understood on 
the unit cell e.g. 9 G [0, 27 r) or on the lift, e.g. 0 G M modulo 27 r. 


5.1 Isochronal coordinates 

Consider 0 with e = 0. The asymptotic (or latent) phase of a point xq in the basin of attraction 
of the limit cycle T of period T is the value of 'd(xo) such that 

lim |x(t) — u{t + '!?(xo))| = 0 , 

t^OO 

where x{t) is a trajectory starting at xq. Thus if u{t) and x{t) are trajectories on and off the 
limit cycle respectively, they have the same asymptotic phase if the distance between u{t) and x{t) 
vanishes as t —)> 00 . The locus of all points with the same asymptotic phase is called an isochron. 
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Symbol 

Phase 

space 

Period 

Uncouplec 

equation 

Advantages 

Disadvantages 


M/TZ, 

[o,r) 

T 


Simplicity of 
uncoupled 
equation, in¬ 
terpretation of 
■Q as “time” 

Phase space 
depends on 

parameters 
and initial 

conditions. 

e 

[0,27r) 

27r 


Phase space 
fixed, good for 
heterogeneous 
oscillators. 

Equations 
need scaling. 

(P 

M/Z, 

[0,1) 

1 


Phase space 
fixed, good for 
heterogeneous 
oscillators. 

Equations 
need scaling. 


Table 2: Comparison of the three conventions for a phase variable that we use in this review. The 
(j) is used for IF models. 


Thus an isochron extends the notion of phase off the cycle (within its basin of attraction). Isochrons 
can also be interpreted as the leaves of the stable manifold of a hyperbolic limit cycle. They fully 
specify the dynamics in the absence of perturbations |127] . 

There are very few instances where the isochrons can be computed in closed form (though see 
the examples in |351j for plane-polar models where the radial variable decouples from the angular 
one). Computing the isochron foliation of the basin of attraction of a limit cycle is a major challenge 
since it requires knowledge of the limit cycle and therefore can only be computed in special cases 
or numerically. 

One computationally efficient method for numerically determining the isochrons is backward 
integration, however it is unstable and in particular for strongly attracting limit cycles the trajecto¬ 
ries determined by backwards integration may quickly diverge to infinity. See Izhikevich |146] for a 
MATLAB code which determines smooth curves approximating isochrons. Other methods include 
the continuation based algorithm introduced by Osinga and Moehlis |243) . the geometric approach 
of Guillamon and Huguet to find high order approximations to isochrons in planar systems |129j . 
quadratic and higher order approximations |281[ 1314) , and the forward integration method using 
the Koopman operator and Fourier averages as introduced by Mauroy and Mezic |212j . This latter 
method is particularly appealing and given its novelty we describe the technique below 

The Koopman operator approach for constructing isochrons for a T-periodic orbit focuses on 
tracking observables (or measures on a state space) rather than the identification of invariant sets. 
The Koopman operator, K, is dehned hy K = z o <I>i(x), where z : M” —)> M is some observable of 
the state space and ‘^^(x) denotes the flow evolved for a time t, staring at a point x. The Fourier 
average of an observable z is defined as 


z{x; Lu) 


lim ^ f {zo^t){x)e 
T—^oo 1 Jq 


(18) 


For a fixed x, (18) is equivalent to a Fourier transform of the (time-varying) observable computed 
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along a trajectory. Hence, for a dynamics with a stable limit cycle (of frequency 2tt jT), it is clear 
that the Fourier average can be nonzero only for the frequencies ojn = n G Z. The Fourier 

averages are the eigenfunctions of K, so that 

K'z{x;ujn) = z{x-,u!n), n G Z. 


Perhaps rather remarkably the isochrons are level sets of z{x;ujn) for almost all observables. The 
only restriction being that the first Fourier coefficient of the Fourier observable evaluated along the 


limit cycle is nonzero over one period. An example of the use of this approach is shown in Fig. 11 
where we plot the isochrons of a Stuart-Landau oscillator. 



Figure 11: Isochrons found using the method of Fourier averages for the Stuart-Landau oscillator: 

= z[\{l + ic)/2 + ioj] — 2;|zp(l + ic)/2, z = x + iy, with A = 2, c = 1 and a; = 1 so that T = 27r. 
The black curve represents the periodic orbit of the system; in this case, we have 6 = "d. The 
background colour represents the Fourier average, whilst the coloured lines are the isochrons, given 
as level sets of the Fourier average. The white dots are the actual isochrons, computed analytically 
as (x, y) = ((1 + r) cos{9 + cln(l -|- r)), (1 + r) sin(0 -|- cln(l -|- r))), r G (—1, oo). 


5.2 Phase-amplitude models 


An alternative (non isochronal) framework for studying oscillators with an attracting limit cycle is 
to make a transformation to a moving orthonormal coordinate system around the limit cycle where 
one coordinate gives the phase on the limit cycle while the other coordinates give a notion of distance 
from the limit cycle. It has long been known in the dynamical systems community how to construct 
such a coordinate transformation, see |133j for a discussion. The importance of considering the 
effects of both the phase and amplitude interactions of neural oscillators has been highlighted 
bv several authors including Ermentrout and Kopell |104j and Medvedev m, and that this is 
especially pertinent when considering phenomenon such as oscillator death (and see § 3.3). Phase- 


amplitude descriptions have already successfully been used to find equations for the evolution of the 
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energies (amplitudes) and phases of weakly coupled weakly dissipative networks of nonlinear planar 
oscillators (modelled by small dissipative perturbations of a Hamiltonian oscillator) [T71 [IHl Il9] . 
Lee et al. M use the notion of phase and amplitudes of large networks of globally coupled 
Stuart-Landau oscillators to investigate the effects of a spread in amplitude growth parameter 
(units oscillating with different amplitudes and some not oscillating at all) and the effect of a 
homogeneous shift in the nonlinear frequency parameter. 

We now discuss the phase-amplitude coordinate transformation detailed by Hale |133j . some of 
the situations where it has been employed and other areas in which a phase-amplitude description 
of oscillators is necessary to reveal the dynamics of the system. Consider again the system 0 
with e = 0 which has an attracting hyperbolic periodic orbit. Make a transformation to a moving 
orthonormal coordinate system as follows. Choose one axis of the coordinate system to be in the 
direction of unit tangent vector along the periodic orbit, given by 




du 


di? 


The remaining coordinate axes can be expressed as the columns of an n x (n — 1) matrix We 
can then write an arbitrary point x in terms of its phase i? G [0, T) and its amplitude p: 

x{'d,p) = u{'&) + C{'&)p. 

Here, \p\ represents the Euclidean distance from the limit cycle. The vector of amplitudes p G 
allows us to consider points away from the limit cycle. 

Upon projecting the dynamics onto the moving orthonormal system, we obtain the dynamics 
of the transformed system for G [0, T) and p G 

= l + /i(i?,p), = A{^)p +f2{^,p). 


where 


h{^,p) = + [f{u + C,p) - fiu )], 

/ 2 (??, p) = ~ •/’('“) “ ^f(P\ > 


h{'d,p) = 


du 

d^ 


I tT dC 


-1 


e, A{^)=e 




(19) 


and D/ is the Jacobian of the vector field /, evaluated along the periodic orbit u. The technical 
details to specify the orthonormal coordinates forming C, can be found in the appendix of [344] . It 
is straight-forward to show that fi{d,p) —> 0 as \p\ —> 0, /2(i?, 0) = 0 and that df 2 {'&, 0 )/dp = 0. 
In the above, ^(1?) describes the -(^-dependent rate of attraction or repulsion from cycle and /i 
captures the shear present in the system, that is, whether the speed of -i? increases or decreases 
dependent on the distance from cycle. A precise definition for shear is given in |245| and will be 
discussed further in 


5.3 


Some caution must be exercised when applying this transformation as it will break down when 
the determinant of the Jacobian of the transformation vanishes. This never occurs on cycle (where 
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p = 0) but it may do so for some \p\ = k > 0, setting an upper bound on how far from the limit 
cycle these phase-amplitude coordinates can be used to describe the system. In |344| it is noted 
that for the planar ML model the value of k can be relatively small for some values of •d, but that 
breakdown occurs where the orbit has high curvature. In higher dimensional systems this issue 
would be less problematic. 

Similarly, the coordinate transformation can be applied to driven systems (i.e. © with e 7^ 0) 
where e is not necessarily small. In this case the dynamics in (-d, p) coordinates, where G [0, T) 
and p G are 

= 1 + p) + eh^{^, p)g{u{^) + C{^)p, t), (20) 

Ap = A{d)p + /2('i9, p) + p)g{u{d) + C(i?)p, t), (21) 

with 

B{d,p)=ln-^ph^{^,p), (22) 

and In is the n x n identity matrix. Here, h and B describe the effect in terms of "d and p that the 
perturbations have. For planar models, H = I2. 

Medvedev m has employed this phase-amplitude description to determine conditions for 
stability of the synchronised state in a network of identical oscillators with separable linear coupling. 
Medvedev |215j has also used the framework to consider the effects of white noise on the synchronous 
state, identifying the types of linear coupling operators lead to synchrony in a network of oscillators 
provided that the strength of the interactions is sufficiently strong. 


5.3 Dynamics of forced oscillators: shear-induced chaos 


Since phase-amplitude coordinates can capture dynamics a finite distance away from the limit cycle, 
(and additionally have the advantage over isochronal coordinates of being defined outside of the 
basin of attraction of the limit cycle) they can be used to model dynamical phenomena in driven 
systems where the perturbations necessarily push the dynamics away from the limit cycle. There 
is no need to make any assumptions about the strength of the forcing e. 

The phase-amplitude description of a forced oscillator is able to detect the presence of other 
structures in the phase space. For example if the system were multistable, phase-amplitude co¬ 
ordinates would track trajectories near these other structures and back again, should another 
perturbation return the dynamics to the basin of attraction of the limit cycle. These coordinates 
would also detect the presence of other non-attracting invariant structures such as saddles in the 
unperturbed flow. Orbits passing near the saddle will remain there for some time and forcing may 
act to move trajectories near to this saddle before returning to the limit cycle. It may also be 
the case that the forcing acts to create trapping regions if the forcing is strong compared to the 
attraction to the limit cycle. 

Another dynamical phenomenon which phase-amplitude coordinates are able to capture is the 
occurrence of shear-induced chaos. Shear in a system near a limit cycle F is the differential in 
components of the velocity tangent to the limit cycle as one moves further from F in phase space. 
Shear forces within the system act to speed up (or slow down) trajectories further away from the 
limit cycle compared to those closer to it. This phenomenon is illustrated in Fig. 12 


As we show below, an impulsive force is applied (the system is kicked) a ‘bump’ in F is produced. 
If there is sufficient shear in the system then the bump is folded and stretched as it is attracted 
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Figure 12: The stretch-and-fold action of a kick followed by relaxation in the presence of shear. 
Left: A non-constant kick moves data away from the limit cycle (horizontal grey line). As the 
image of the cycle under the kick relaxes back to the cycle (under the unpertnrbed flow) the action 
of shear causes folds to appear. Right: The geometry of folding in relation to the isochron foliation 
(black lines). An initial segment 70 is kicked to the blue curve and is then allowed to relax back 
to cycle, passing throngh the red and green curves, which are images of the blue curve under the 
nnperturbed flow. 


back to the limit cycle. Such folding can potentially lead to the formation of horseshoes and strange 
attractors. However, if the attraction to the limit cycle is large compared to the shear strength or 
size of the kick then the bumps will dissipate before any significant stretching occurs. 

Shear-induced chaos is most commonly discussed in the context of discrete time kicking of limit 
cycles. Wang and Young |339l 134011338j prove rigorous results in the case of periodically kicked 
limit cycles with long relaxation times. Their results provide details of the geometric mechanism 
for producing chaos. Here we briefly review some of these results. More detailed summaries can be 
found in m and m\- 

Let be the flow of the unpertnrbed system, which has an attracting hyperbolic limit cycle 
T. We can think of a kick as a mapping k. The dynamics of the system with periodic kicking 
every r units of time can be obtained by iterating the map TV = <!>,- o k . Provided that there is a 
neighbourhood ZV of T such that points in U remain in the basin of attraction of T under kicks and 
r is long enough that the kicked points can return to U, then T,- = nn>o-FV(ZV) is an attractor for 
the periodically kicked system. If the kicks are weak then T,- is generally expected to be a slightly 
perturbed version of T and we should expect fairly regular behaviour. In this case Ft- is an invariant 
circle. To obtain more complicated dynamics it is necessary to break the invariant circle. This idea 
is illustrated by the following linear shear model considered in [MO]: 

l« = i + ’yy. 

A ^ 

—y = -Ay -k AH{'d) ^ 6{t + nr), 

n=0 

where (??, y) G [0,T) x R = x R are coordinates in the phase space, X, a, A > 0 are constants 
and H —)> M is a nonconstant smooth function. The unforced system {A = 0) has a limit cycle 

F = X {0}. If the quantity 

<7 , shear 

—A =- X (kick amplitude ), 

A contraction rate 

is sufficiently large, then it is possible to show that there is a positive measure set T C such 
that for all r G T, F,- is a strange attractor of Fr |338j . How large the quantity must be depends 
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on the function Since H(-d) is nonconstant (to create the bumps), the larger shear a and 

the kick amplitude A, the more folding will occur as the bump is attracted back to the limit cycle. 
Also note that weaker limit cycles (smaller A) create more favourable conditions for shear induced 
chaos as the shear has longer to act before the bump is attracted back to the limit cycle. 

For a more general system with periodic forcing the shear may not appear explicitly as a 
parameter. To elucidate what kind of kicks may cause shear induced chaos in this case we appeal 
to the isochrons of the system. We think of the isochrons of F as the strong stable manifold of F. 
Suppose, as illustrated in Fig. [T^ that a section Fq of the limit cycle is kicked upwards with the 


end points held fixed and assume t = np for some n,p € TLA. Since the isochrons are flow invariant, 
during relaxation the flow moves each point of the kicked curve k{Tq) back towards F along the 
isochrons. In Fig. 12 we can clearly see the effect of the shear with a fold forming. 

From Fig. 12 one can see that kicks along isochrons or in directions roughly parallel to the 


isochrons will not produce strange attractors, nor will kicks that carry one isochron’s manifold 
to another. The cause of the stretching and folding is the variation in how far points x G F are 
moved by k in the direction transverse to the isochrons (i.e. the ordering of points in terms of 
asymptotic phase is altered by the action of the kick). Lin and Young |196j emphasise that the 
occurrence of shear induced chaos depends on the interplay between the geometries of the kicks 
and the dynamical structures of the unforced system. 

In the case of the linear shear model above, the isochrons of the unforced system are simply 
the lines with slope —Xja in coordinates. Variation in kick distances in directions transverse 

to these isochrons is guaranteed with any nonconstant function H, with greater variation given by 
larger values of cr/A and A. 

Beyond the rigorous results proved by Wang and Young |339l I340L133811196j concerning peri¬ 
odically kicked limit cycles of the linear shear model and for supercritical Hopf bifurcations, Ott 
and Stenlund |245j prove that shear-induced chaos may exist near general limit cycles. In addition, 
Lin and Young |196) have carried out numerical studies including random kicks at times given by a 
Poisson process and systems driven by white noise. They also consider forcing of a pair of coupled 
oscillators. In all cases, shear-induced chaos occurs when the shearing and amplitude of the forcing 
are large enough to overcome the effects of damping. 

Lin et al. dM! demonstrate that the ML model can exhibit shear-induced chaos near the 
homoclinic bifurcation when periodically forced, by plotting images of the periodic orbit under 
successive applications of the kick map and calculating the maximum Lyapunov exponent. They 
also emphasise that the phenomenon of shear-induced chaos cannot be detected by the perturbative 
techniques such as those outlined in 


5.4 and § 5.5 below. 


Wedgwood et al. [344] go on to show that the phase-amplitude description is well-suited to 
understanding the behaviour of neural oscillators in response to external stimuli. They consider a 
number of neural oscillator models in addition to a generic planar model and examine the response 
of the system to periodic pulsatile forcing by taking x G M^, 


9{x, t) = - nr), 0)^, 

and a variety of choices for /(x) in ©• Their numerical simulations indicate that for both 
the FHN and ML models the behaviour remains qualitatively similar when the relevant functions 
are replaced by ap, f 2 {'&,p) is dropped and A{'d) is replaced with —A (for a wide range of 
cr, A > 0). They then focus on the effect of the form of the forcing function in the phase-amplitude 
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coordinates. Evaluating the functions and of equations (19) and (22) respectively 


for the particular model and denoting the first component of h as Pi and the first component of ^ as 
P 2 they find that forcing each system at the same ratio of the natural frequency of the underlying 
periodic orbit and implementing the choices above gives the following ODE on G [0, T), p G 


dt 


'd = I + ap + ePi{'d,p) ^ - nr), 


Ap = -Xp + eP 2 (i?) ^ - nr). 


By developing a stroboscopic map for this system, Wedgwood et al. [344j numerically evaluate the 
largest Lyapunov exponent which is found to be positive for the Morris-Lecar model but negative 
for the FHN model. 

The analysis of the behaviour of generic networks of oscillators within a phase-amplitude frame¬ 
work is a challenging open problem but such a description would allow for greater accuracy (com¬ 
pared to the phase only methods traditionally used and described below) in elucidating a richer 
variety of the complex dynamics of oscillator networks. 


5.4 Phase oscillator models 


To obtain a phase description, one can consider the limit of strong attraction in equations (20)-(21). 


However, it is more appealing to consider a phase variable with a uniform rotation rate that assigns 
a phase coordinate i? G [0, T) to each point a: G T according to 


dt 


'&{x{t)) = 1, 


X G r. 


(23) 


This reduction to a phase description gives a simple dynamical system, albeit one that cannot 
describe evolution of trajectories in phase-space that are away from the limit cycle. However, the 
phase reduction formalism is useful in quantifying how a system (on or close to a cycle) responds 
to weak forcing, via the construction of the inhnitesimal phase response curve (iPRC). This can 
be written, for a given ODE model, as the solution to an adjoint equation. For a given high 
dimensional conductance based model this can be solved for numerically, though for some normal 
form descriptions closed form solutions are also known [55]. 

The iPRC at a point on cycle is equal to the gradient of the (isochronal) phase at that point. 
Writing this vector quantity as Q means that weak forcing eg{t) in the original high dimensional 
models transforms as did/dt = 1 -|- e{Q,g) where (•,•) defines the standard inner product and e 
is a small parameter. For periodic forcing such equations can be readily analysed, and questions 
relating to synchronisation, mode-locking and Arnol’d tongues can be thoroughly explored |255j . 
Moreover, this approach forms the basis for constructing models of weakly interacting oscillators, 
where the external forcing is pictured as a function of the phase of a firing neuron. This has led to 
a great deal of work on phase-locking and central pattern generation in neural circuitry and see for 
example |142| . 

However, the assumption that phase alone is enough to capture the essentials of neural response 
is one made more for mathematical convenience than being physiologically motivated. Indeed for 
the popnlar type I ML firing model with standard parameters, direct numerical simulations with 
pulsatile forcing show responses that cannot be explained solely with a phase model [194] . as just 
highlighted in § 5.3 (since strong interactions will necessarily take one away from the neighbourhood 


of a cycle where a phase description is expected to hold). 
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5.5 Phase response 

It is common practice in neuroscience to characterise a neuronal oscillator in terms of its phase 
response to a perturbation [573]. This gives rise to the notion of a so-called phase response curve 
(PRC), which for a real neuron can be determined experimentally |114[I317] . and can also be related 
to the poststimulus time histogram |130j . Consider a dynamical system dx/dt = f{x), x G 
with a T-periodic solution u{t) = u{t + T) and introduce an infinitesimal perturbation Axq to the 
trajectory u{t) at time t = 0. This perturbation evolves according to the linearised equation of 
motion: 

—Ax = Df{u{t))Ax, Ax(0) = Axq. 

Here Df{u) denotes the Jacobian of / evaluated along u. Introducing a time-independent isochronal 
phase shift At} as + Au{t)) — we have to first order in Ax that 

Ai? = (Q(t),Ax(t)), (24) 


where (•, •) defines the standard inner product, and Q = is the gradient of iD evaluated at u{t). 


Taking the time-derivative of (24) gives 


Q, Ax^ = - (^Q, = - {Q, Df{u)Ax) = - {Df'^{u)Q, Ax). 


Since the above equation must hold for arbitrary perturbations, we see that the gradient Q = Vu'& 
satishes the linear equation 

= D{t)Q, D{t) = -Df{u{t)), (25) 

subject to the boundary conditions 


(V^.(o)^?,/(M(0))) = 1 and Q{t) = Q{t + T). 


The first condition simply guarantees that d'd/dt = 1 (at any point on the periodic orbit), and 
the second enforces periodicity. Note that the notion of phase response can also be extended to 
time-delayed systems |168( I236| . 


In general equation (25) must be solved numerically to obtain the iPRC, say, using the adjoint 
routine in XPPAUT |99| or MatCont |126| . However, for the case of a nonlinear IF model, defined 
by (§, the PRC is given explicitly by 


Qi^) = 


1 


/ + /oT-i(79)’ 


pvt 

T(x) = r / 

■tvB. 


HP du' 


I + f{v') 


A ’ 


(26) 


where the period of oscillation is found by solving 'l'(xth) = T (and the response function is valid for 
finite perturbations). For example, for the quadratic IF model, obtained from ([^ with the choice 
with f{v) = v^, and taking the limit Xth —^ oo and xr —>■ —oo we hnd Q('i?) = sin^('d7r/r)/I with 
T = TTr/y/l, recovering the iPRC expected of an oscillator close to a SNIC bifurcation [Ml 1107] . 
For a comprehensive discussion of iPRCs for various common oscillator models see the excellent 
survey article by Brown et al. [55]. The iPRC for planar piecewise linear IF models can also be 


computed explicitly although we do not discuss this further here. In Fig. we show the 
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Figure 13: From left to right: Phase response curves for the Hodgkin-Huxley model, the Wilson 
model and the Morris-Lecar model corresponding to the discussion in § [^ Orbits are in blue and 
iPRCs in red. 


numerically computed iPRCs for the Hodgkin-Huxley model, the Wilson cortical model and the 
Morris-Lecar model previously discussed in § [^ For completeness it is well to note the iPRC for 
a planar oscillator close to a supercritical Hopf bifurcation has an adjoint with a shape given by 
(—sin(27rt/r), cos(27rt/r)) for an orbit shape (cos(27rf/r), sin(27rt/r)). 

Having dehned the phase in some neighbourhood of the limit cycle, we can write (23) as 


dt 


'd{x) 



For dx/dt = f{x) this gives {'Vx'&, f{x)) = 1. We now consider the effect of a small external periodic 
force on the self sustained oscillations as in with g{x,t) = g{x,t + A), where in general A is 
different to T. For weak perturbations (e <C 1) the state point will leave F, though will stay in some 
neighbourhood, which we denote by U. We extend the phase off cycle using isochronal coordinates 
so that {Vx'&, f{x)) = 1 holds for any point x gU. For a perturbed oscillator, in these coordinates 
we have 

^x'^ = {Vx'&, fix) + egix, t)) = l + e g{x, t)). 

As a first approximation we evaluate the right hand side of the above on the limit-cycle to get an 
equation for the phase dynamics: 


—d = l + eli'd,t), 


Ii'&,t) = {Qiu{i^)),giui^),t)). 


(27) 


The phase dynamics (27) is still very hard to analyse, and as such it is common to resort to a 


further simplification, namely averaging. First let us introduce a rotating phase ^j; = i!} — Tt/A, in 


which case (27) becomes 


dt 


V’ 


-6 + el{ip -I- Tt/A,t), 


6 = T/A-l. 


If both e and 6 are small then d'ljj/dt ~ 0 and evolves slowly. Thus we may set “ipis) ~ 'ipit) for 
s G [t,t + T]. Averaging the above over one period gives 


d 

dt 


~ — d -|- eHitp), 



{Qiui^p + s)),giui'if + s), s)) ds. 


(28) 
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where we have used the result that I{'4) + s,t) = /(V' + s + T, t+T). The function is T-periodic 
and can be written as a Fourier series , with the simplest example of an 

averaged phase-dynamics being 

—ijj = —5 + esin'0 

at 

which is called the Adler equation [6]. If we denote the maximum and minimum of H by -fimin 
and -fimax respectively then for a phase-locked 1:1 state defined by = 0 we require eH^\^ < 
6 < ei/max- In this case there are two fixed points ’ip± defined by = 6. One of these is 

unstable (say V'-, so that eH'{'tp_) > 0) and the other stable (■0+, with ei7'('0+) < 0). This gives 
rise to a rotating solution with constant rotation frequency so that i? = (1 + 5)t + The two 
solutions coalesce in a saddle-node bifurcation when 5 = eflmin and 5 = e ffmax (or equivalently 
when H'{'ip±) = 0). In the case of the Adler model the parameter region for phase-locking is given 
explicitly by a triangular wedge defined by e = |5| - a so-called Arnol’d tongue. Outside of this 
tongue solutions drift (they cannot lock to the forcing period) according to ( [28| ), and the system 
evolves quasi-periodically. We treat weakly coupled phase oscillators in § [^ 


6 Weakly coupled phase oscillator networks 


The theory of weakly coupled oscillators |173( 1102] is now a standard tool of dynamical systems 
theory and has been used by many authors to study oscillatory neural networks, see for example 
[MIIMIIIMIM. The book by Hoppensteadt and Izhikevich provides a very comprehensive review 
of this framework [142] . which can also be adapted to study networks of relaxation oscillators (in 
some singular limit) [1471171] . 

Consider, for illustration, a system of interacting limit-cycle oscillators (|8|). Following the 
method in 


5.5, similar to (28) we obtain the network’s phase dynamics in the form 


dt 


N 


9i = LOi + e'^Wij{Qi{ui{6i)),G{uj{9j))), 
i=i 


(29) 


where the frequency Wj allows for the fact that oscillators are not identical and, for this reason, 
we will assume that 9i G [0,27r). Precisely this form of network model was originally suggested 
by Winfree to describe populations of coupled oscillators. The Winfree model [352] assumes a 
separation of time-scales so that an oscillator can be solely characterised by its phase on cycle (fast 
attraction to cycle) and is described by the network equations 


1 1 _ 

-9,=u;i + eR{9i)-Y,m)^ ( 30 ) 

i=i 

describing a globally coupled network with a biologically realistic PRC R and pulsatile interaction 
function P. Using a mixture of analysis and numerics Winfree found that with large N there was 
a transition to macroscopic synchrony at a critical value of the homogeneity of the population. 
Following this Kuramoto [173] introduced a simpler model with interactions mediated by phase 
differences, and showed how the transition to collective synchronisation could be understood from 
a more mathematical perspective. For an excellent review of the Kuramoto model see [312] and [5]. 
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The natural way to obtain a phase-difference model from (29) is, as in § 5.5, to average over one 


period of oscillation. For simplicity let us assume that all the oscillators are identical, and uji = u 
Vh in which case we find that 


dt 


9i=u} 


N 

+ e’^WijHiOj - 9i), 
i=i 


(31) 


where 


= {Qiu{s)),G{u{'ip + s)))ds. 

The 27r-periodic function H is referred to as the phase interaction function. 
Fourier series for Q and G as 


(32) 


If we write complex 


g(t) = ^g„e‘"* and G(t) = 


respectively then 


HW = '^Hne 


in'll) 


(33) 




with Hn = {Q-n, Gn). Note that certain caution has to be exercised in applying averaging theory. 
In general, one can only establish that a solution of the unaveraged equations is e-close to a 
corresponding solution of the averaged system for times of No such problem arises in the 

case of hyperbolic fixed points corresponding to phase-locked states. 

When describing a piece of cortex or a central pattern generator circuit with a set of oscillators, 
the biological realism of the model typically resides in the phase interaction function. The simplest 
example is H(ip) = sin('0), which when combined with a choice of global coupling defines the well 
known Kuramoto model mi- However, to model realistic neural networks one should calculate 
(32) directly, using knowledge of the single neuron iPRC and the form of interaction. As an 

that can be written in the form G{u{'if)) = 


2.5 


example consider synaptic coupling, described in 

+ W'2vr), and a single neuron model for which the iPRC in the voltage variable v is given 
by R (say experimentally or from numerical investigation). In this case 


= / R{2'ks — V’)?7(27rs)ds. 

Jo 


(34) 


If instead we were interested in diffusive (gap-junction) coupling then we would have 

H{'P) = — / R{s)[v{s+ ip) - v{s)]ds. 

27r Jo 

For the HH model R{t) is known to have a shape like — sin(t) for a spike centred on the origin (see 
Fig. 13). Making the further choice that r/(t) = a^te~°‘^ then (34) can be evaluated as 




[1 — (1/a)^] sin('0) — 2/acos(V') 

27r[l + (l/a)2]2 • 


(35) 


43 










In the particular case of two oscillators with reciprocal coupling and w = 1 then 


^0i = l + e/7(02-0i), 
^02 = l + e//(0i-02), 


and we define (p := 02(i) ~ Qi{t)- A phase-locked solution satisfies constant phase difference (f) that 
is a zero of the (odd) function 

K{p) = e[H{-p)-H{ip)]. 


A given phase-locked state is then stable provided that K'{p) < 0. Note that by symmetry both 
the in-phase {p = 0) and anti-phase {p = tt) states are guaranteed to exist. For the form of phase 
interaction function given by (35), The stability of the synchronous solution is governed by the sign 
of i^'(O): 

sgn iF'(O) = sgn {-eF'(O)} = sgn {-e[(l - (1/a))^]} . 


Thus for inhibitory coupling (e < 0) synchronisation will occur if l/a > 1, namely when the synapse 
is slow (a —>■ 0). It is also a simple matter to show that the anti-synchronous solution {p = vr) is 
stable for a sufficiently fast synapse {a oo). It is also possible to develop a general theory for 
the existence and stability of phase-locked states in larger networks than just a pair. 


6.1 Phase, frequency and mode locking 

Now suppose we have a general population of N > 2 coupled phase oscillators 

described by phases 9j G M/27rZ. For a particular continuous choice of phases 6{t) for the trajectory 
one can define the frequency of the jth oscillator as 

% = f[^j{T) - 0j(O)]. 

This limit will converge under fairly weak assumptions on the dynamics |154| . though it may vary 
for different attractors in the same system, for different oscillators j and in some cases it may vary 
even for different trajectories within the same attractor. 

We say two oscillators j and k are phase locked with ratio (n : m) for (n, m) £ 7? \ (0,0) with 
no common factors of n and m, if there is an M such that 

\n9j{t) — m9k{t)\ < M, 

for all t > 0. The oscillators are frequency locked with ratio (n : m) if 

nVtj — mQk = 0 - 

If we say they are simply phase (or frequency locked) without explicit mention of the {n : m) ratio, 
we are using the convention that they are (I : 1) phase (or frequency) locked. The definition of 
Ilj means that if two oscillators are phase locked then they are frequency locked. The converse is 
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not necessarily the case: two oscillators may be frequency locked but not phase locked if the phase 
difference nBj{t) — m6k{t) grows sublinearly with t. 

For the special case globally coupled averaged networks {wij = 1/N for the system ([^) is S^xT 
equivariant. By topological arguments, maximally symmetric solutions describing synchronous, 
splay, and a variety of cluster states are expected to exist generically [23]. The system © with 
global coupling is in itself an interesting subject of study in that it is of arbitrarily high dimension 
N but is effectively determined by the single function H{(p) that is computable from a single pair 
of oscillators. The system (and variants thereof) have been productively studied by thousands of 
papers since the seminal work of Kuramoto [173j . 


6.2 Dynamics of general networks of identical phase oscillators 


The collective dynamics of phase oscillators have been investigated for a range of regular network 
structures including linear arrays and rings with uni- or bi-directional coupling e.g. [102119711231IM] . 
and hierarchical networks |295j . In some cases the systems can be usefully investigated in terms 


of permutation symmetries of (31) with global coupling, for example Z^r or for uni- or bi¬ 
directionally coupled rings. In other cases a variety of approaches have been developed adapted 
to particular structures though these have not in all cases been specifically applied to oscillators; 
some of these approaches are discussed in 


3.2 


We recall that the form of the coupling in (31) is special in the sense that it assumes the 
interactions between two oscillators are independent of any third - pairwise coupling 
there are degeneracies such as 


ED. If 


m—1 

E 

fc =0 




2'Rk 


m 


= 0 , 


(36) 


which can appear when some of the Fourier components of H are zero, this can lead to degeneracies 
in the dynamics. For example |342| . while |23[ Theorem 7.1] shows that if H satisfies (36) for 


some m > 2 and is a multiple of m then (31), with global coupling, will have m-dimensional 


invariant tori in phase space that are foliated by neutrally stable periodic orbits. This degeneracy 
will disappear on including either non-pairwise coupling or introducing small but non-zero Fourier 
components in the expansion of H but as noted in |167j this will typically be the case for the 
interaction of oscillators even if they are near a Hopf bifurcation. 

We examine in more detail some of the phase locked states that can arise in weakly coupled 


networks of identical phase oscillators described by (31). We define a 1:1 phase-locked solution to 


be of the form 9i{t) = ipi + Qt, where (fi is a constant phase and 11 is the collective frequency of 


the coupled oscillators. Substitution into the averaged system (31) gives 


N 


n = w -I- e © WijH{ipj - ifi). 


(37) 




After choosing some reference oscillator, these N equations determine the collective frequency H 
and N — 1 relative phases with the latter independent of e. 

It is interesting to compare the weak coupling theory for phase-locked states with the analysis 
of LIF networks from § 4.3 Equation (13) has an identical structure to that of equation (37) 


(for li = I for all i), so that the classification of solutions using group theoretic methods is the 
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same in both situations. There are, however, a number of significant differences between phase¬ 
locking equations (37) and (13). First, equation (13) is exact, whereas equation (37) is valid only 
to 0(e) since it is derived under the assumption of weak coupling. Second, the collective period of 


oscillations A must be determined self-consistently in equation (13). 


In order to analyse the local stability of a phase-locked solution $ = (c;/!)!,..., we linearise 
the system by setting 9i{t) = ipi + Vtt + 9i{t) and expand to first-order in 0*: 


N 


1=1 


where 


N 

= WijH'{ipj - cpi) - 6ij'^WikH'{ipk - (fi), 

k=l 

and = dH{ip)/dip. One of the eigenvalues of the Jacobian T-L is always zero, and the cor¬ 

responding eigenvector points in the direction of the flow, that is (1,1,..., 1). The phase-locked 
solution will be stable provided that all other eigenvalues have a negative real part. We note that 
the Jacobian has the form of a graph-Laplacian mixing both anatomy and dynamics, namely it is 
the graph-Laplacian of the matrix with components —WijH'{ipj — p>i). 


6.2.1 Synchrony 

Synchrony (more precisely, exact phase synchrony) is where 9i = 92 = ■ ■ ■ = 9 n-i = 9n = + to 

for some fixed frequency U is a classic example of a phase-locked state. Substitution into (31), 
describing a network of identical oscillators, shows that 11 must satisfy the condition 


N 

n = cj -|- eH (0) w, 

i=i 


ij, 


Mi. 


One way for this to be true for all i is if 77(0) = 0, which is the case say for H[9) = sin(0) or for 
diffusive coupling, which is linear in the difference between two state variables so that 77(0) = 0. 
The existence of synchronous solutions is also guaranteed if '^ij independent of i. This 

would be the case for global coupling where Wij = 1/N, so that the system has permutation 
symmetry. 

If the synchronous solution exists then the Jacobian is given by e77'(0)£ where L is the graph- 
Laplacian with components Cij = Wij — 6ij Wik- We note that £ has one zero eigenvalue, with 
eigenvector (1,1,...,1,1). Hence if all the other eigenvalues of £ lie on one side of the imaginary axis 


then stability is solely determined by the sign of ei7'(0). In Fig. 14 we show the eigenvalues of the 
graph-Laplacian of the anatomical network structure of the Macaque monkey brain, as determined 
from the CoCoMac database [305] . Here all the eigenvalues lie to the left of the imaginary axis so 
that the stability of the synchronous solution (should it exist) is solely determined by the sign of 
e77'(0). For global coupling we have that Cij = N~^ — dij, and the {N — 1 degenerate) eigenvalue 
is —1. Hence the synchronous solution will be stable provided A = —e77'(0) < 0. 
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Figure 14: Eigenvalues of the graph-Laplacian for the anatomical network structure of the Macaque 
monkey brain — generated from the CoCoMac database. 
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6.2.2 Asynchrony 

Another example of a phase-locked state is the purely asynchronous solution whereby all phases 
are uniformly distributed around the unit circle. This is sometimes referred to as a splay state or 
splay-phase state and can be written dOi/dt = with 0j_|_i — 6i = 2 tt/N Mi. Like the synchronous 
solution it will be present but not necessarily stable in networks with global coupling, with an 
emergent frequency that depends on H\ 



In this case the Jacobian takes the form 

^ [^n—m hnmT'] j 

where T = LJ'(27r/c/A) and An = H'{2'Kn/N). Hence the eigenvalues are given by \p = 

e\vp — r]/A, p = 0,..., A — 1 where Up are the eigenvalues of An-m- Ylm ^n-mO?n = i^pO-n, where 
an denote the components of the pth eigenvector. This has solutions of the form an = so 

that Vp = giving 


m=l 


\ rg27rim; 


V A j 


1 ], 


and the splay state will be stable if Re (Ap) < 0 Vp 7 ^ 0. 

In the large A limit A —I- 00 we have the useful result that (for global coupling) network 
averages may be replaced by time averages: 


lim — 
N^oo A 



F{t)dt = To, 


for some 27r-periodic function F{t) = F{t 2 tt) (which can be established using a simple Riemann 
sum argument), with a Fourier series F{t) = Hei^ce in the large A limit the collective 
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frequency of a splay state (global coupling) is given by 0 = w + eH ^, with eigenvalues 


\p = — / = —2iTipeH-r 

271 - Jo 


Hence a splay state is stable if —pelmi^p < 0, where we have used the fact that since H(9) is real, 
then ImH-p = —ImHp. As an example consider the case H{9) = 9(Tr — 9){9 — 2 tt) for 9 G [0, 27r) 


and e > 0 (where H is periodically extended outside [0,27r)). 
the Fourier coefficients (33) as Hn = 6i/n^, so that —pelmHj 


It is straight forward to calculate 
= —6ejp‘^ < 0 Vp / 0. Hence the 
asynchronous state is stable. If we flip any one of the coefficients Hm —Hm then Re Am > 0 and 
the splay state will develop an instability to an eigenmode that will initially destabilise the system 
in favour of an m-cluster. 



Figure 15: (a) A phase interaction given by eH{9) = 9 {'k — 9){9 — 2tt) for 9 G [0, 27r) with complex 
Fourier series coefficients given by Hn = dijr?. The remaining panels show the effect of flipping 
the sign of just one of the co-efficients, namely Hm —>■ —Hm- (b) m = 1 , and the asynchronous 
solution will destabilise in favour of the synchronous solution since H'{Q) >0. (c) m = 2, and the 
asynchronous solution will destabilise in favour of a 2-cluster, (d) m = 3, and the asynchronous 
solution will destabilise in favour of a 3-cluster. Note that only small changes in the shape of H, 
as seen in panels (c)-(d), can lead to a large change in the emergent network dynamics. 


6.2.3 Clusters for globally coupled phase oscillators 

For reviews of the stability of cluster states (in which subsets of the oscillator population synchro¬ 
nise, with oscillators belonging to different clusters behaving differently) we refer the reader to 
[23l[Il9l[5i]; here we use the notation of [242j . If a group of N oscillators is neither fully synchro¬ 
nised nor desynchronised it may be clustered. We say A = {Ai, ..., Am} is an M-cluster partition 
of{l,...,iV} if 

M 

{l,...,iV}= U Ap, 

p=i 

where Ap are pairwise disjoint sets {Ap n Ag = 0 if p 7 ^ g). NB if Up = \ Ap\ then 

M 

Up=A7. 

p=i 

One can refer to this as a cluster of type (ai,..., um)- It is possible to show that any clustering 
can be realised as a stable periodic orbit of the globally coupled phase oscillator system |242| for 
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suitable choice of phase interaction function; more precisely, there is a coupling function H for the 


system (31), with global coupling, such that for any N and any given M-cluster partition A of 


{1,..., A^} there is a linearly stable periodic orbit realising that partition (and all permutations of 
it). See also [294] where the authors consider clustering in this system where H{(p) = sin Mip. More 
generally, there are very many invariant subspaces corresponding to spatio-temporal clustering that 
we can characterise in the following Theorem. 


Theorem 6.1 [23], Theorem 3.1] The subsets of that are invariant for (31), with global 
coupling, because of symmetries o/Sat x T correspond to isotropy subgroups in the conjugacy class 
of 

‘— {^ki X • • • X Xg '^rn 

where N = mk, k = ki + ■ ■ ■ + ki and Xg denotes the semidirect product. The points with this 
isotropy have im clusters that are split into i groups of m clusters of the size ki. The clusters 
within these groups are cyclically permuted by a phase shift of 27r/m. The number of isotropy 
subgroups in this conjugacy class is ... k^l)]. 

It is a nontrivial problem to discover which of these subspaces contain periodic solutions. Note 
that in-phase corresponds to I = m = 1, ki = N while splay phase corresponds to £ = /ci = 1, 
m = N. The stability of several classes of these solutions can be computed in terms of properties 
of H[(p); see for example 


6.2.1 and § 6.2.2 and for other classes of solution [2311541 [242] . 


6.2.4 Generic loss of synchrony in globally coupled identical phase oscillators 


Bifurcation properties of the globally coupled oscillator system (31) on a varying parameter that 
affects the coupling H{(p) are surprisingly complicated because of the symmetries present in the 
system; see § 3.6 In particular, the high multiplicity of the eigenvalues for loss of stability of the 


synchronous solution means: 


• Path following numerical bifurcation programmes such as AUTO or XPPAUT need to be 
done with great care when applying to problems with N > 3 identical oscillators - these 
typically will not be able to find all solutions branching from one that loses stability. 

• A large number of branches with a range of symmetries may generically be involved in the 
bifurcation; indeed, all 2-cluster states Sk x S^-k- 

• Local bifurcations may have global bifurcation consequences owing to the presence of connec¬ 
tions that are facilitated by the nontrivial topology of the torus [30l |23] . 

• Branches of degenerate attractors such as heteroclinic attractors may appear at such bifur¬ 
cations for > 4 oscillators. 


Hansel et al. [135] consider the system (31) with global coupling and phase interaction function 
of the form 

= — sin((^-|-a)-|-r sin(2(/p), (38) 


for (r, a) fixed parameters; detailed bifurcation scenarios in the cases A = 3,4 are shown in m- 
As an example, Figure shows regions on stability of synchrony, splay phase solutions and robust 
heteroclinic attractors as discussed later in § [^ 
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Figure 16: Regions of stability for the globally coupled N = 4-oscillator system (31) with phase 
interaction function (38) and parameters in the region r>0,0<a<7r m- The narrow stripes 
show the region of stability of synchrony, while the wide stripes show the region of stability of the 
splay phase solution. The pink shaded area shows a region of existence of a robust heteroclinic 
network that is an attractor with in the checkerboard region; the boundaries are described in 
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6.3 Phase waves 


The phase reduction method has been applied to a number of important biological systems, includ¬ 
ing the study of travelling waves in chains of weakly coupled oscillators that model processes such 
as the generation and control of rhythmic activity in central pattern generators (CPGs) underlying 
locomotion |66l I162j and peristalsis in vascular and intestinal smooth muscle |102j . Related phase 
models have been motivated by the observation that synchronisation and waves of excitation can 
occur during sensory processing in the cortex [lOlj . In the former case the focus has been on 
dynamics on a lattice and in the latter continuum models have been preferred. We now present 
examples of both these types of model, focusing on phase wave solutions mi- 


Phase waves: a lattice model 

The lamprey is an eel-like vertebrate which swims by generating travelling waves of neural activity 
that pass down its spinal cord. The spinal cord contains about 100 segments, each of which is a 
simple half-center neural circuit capable of generating alternating contraction and relaxation of the 
body muscles on either side of body during swimming. In a seminal series of papers, Ermentrout 
and Kopell carried out a detailed study of the dynamics of a chain of weakly coupled limit cycle 
oscillators [TMirMirns] . motivated by the known physiology of the lamprey spinal cord. They 
considered N phase-oscillators arranged on a chain with nearest-neighbour anisotropic interactions 
and identified a travelling wave as a phase-locked state with a constant phase difference between 
adjacent segments. The intersegmental phase differences are defined as ipi = 9i+i — Oi. If (/jj < 0 
then the wave travels from head to tail whilst for (pi > 0 the wave travels from the tail to the head. 
For a chain we set Wij = to obtain 



Figure 17: A chain of N phase oscillators pi = 9i+i — 9i with H±{ip) = W±H{ip). 


^^9i=u:i + W+H{92-9i), 

^^9i = iOi + W+H{9i+i - 9i) + W.H{9i-i - 9i), i = 2,...,N-l, 

—= uj^ + W-H{9 n-i - 9n), 

where 9i G Pairwise subtraction and substitution of pi = — 9i leads to an iV — 1 

dimensional system for the phase differences 

= Acui + W+[H{ipi+^) - + W4H{-ipi) - 

for i = 1... A — 1, with boundary conditions H{—ipo) = 0 = H{(pN+i), where Acj, = Ui+i — uji. 
There are at least two different mechanisms that can generate travelling wave solutions. 
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The first is based on the presence of a gradient of frequencies along the chain, that is, Awj 
has the same sign for all i, with the wave propagating from the high frequency region to the low 
frequency region. This can be established explicitly in the case of an isotropic, odd interaction 
function, W± = 1 and H[ip) = —H{—{p), where we have 

= Auji + H{ipi+i) + - 2H{ipi). 

at 

The fixed points $ = {tpi,... ,(Pn) satisfy the matrix equation H{^) = —A~^D, where = 

..., , D = (Awi,..., Alon)'^, and ^ is a tridiagonal matrix with elements An = 

—2, Ai^iJ^i = ^ 1 + 1 ,i = 1. For the sake of illustration suppose that H{ip) = sin((/7 + a). Then a 

solution will exist if every component of A~^D lies between ±1. Let qq = max{|(^“^T))j|}. If 

ao < 1 then for each i = 1,... ,N there are two distinct solutions in the interval [0, 27r] with 
> 0 and II'((p/) < 0. In other words, there are 2^ phase-locked solutions. Linearising 
about a phase-locked solution and exploiting the structure of the matrix A, it can be proven that 
only the solution = {(fi ,..., is stable. Assuming that the frequency gradient is monotonic, 
this solution corresponds to a stable travelling wave. When the gradient becomes too steep to 
allow phase-locking (i.e. ao > 1), two or more clusters of oscillators (frequency plateaus) tend to 
form that oscillate at different frequencies. Waves produced by a frequency gradient do not have a 
constant speed or, equivalently, constant phase lags along the chain. 

Constant speed waves in a chain of identical oscillators can be generated by considering phase- 
locked solutions defined hy (pi = ip for all i, with a collective period of oscillation determined using 
d6*i/dt = n to give Q = LJi + W+H{ipi). The steady state equations are then AuJi + W+H{—ip) = 0, 
Alon-i — W-H{ip) = 0 and Awj = 0, for i = 2,..., A — 2. Thus, a travelling wave solution requires 
that all frequencies must be the same except at the ends of the chain. One travelling solution is 
given by Aww-i = 0 with Awi = —W-H{—ip) and H{ip) = 0. For the choice H{ip) = sin(</5 -|- cr) 
we have that (p = —a and Acji = —IT_ sin(2(T). If 2(T < vr then Aui = UJ 2 — uji < 0 and ui must 
be larger than a ;2 and hence all the remaining uJi for a forward travelling wave to exist. Backward 
swimming can be generated by setting wi = 0 and solving in a similar fashion. 

Phase waves: a continuum model 

There is solid experimental evidence for electrical waves in awake and aroused vertebrate prepa¬ 
rations, as well as semi-intact and active invertebrate preparations, as nicely described by Ermen- 
trout and Kleinfeld m- Moreover, these authors argue convincingly for the use of phase models 
in understanding waves seen in cortex and speculate that they may serve to label simultaneously 
perceived features in the stimulus stream with a unique phase. More recently it has been found 
that cortical oscillations can propagate as travelling waves across the surface of the motor cortex of 
monkeys (Macaca mulatta) |277j . Given that to a first approximation the cortex is often viewed as 
being built from a dense reciprocally interconnected network of corticocortical axonal pathways, of 
which there are typically 10^° in a human brain it is natural to develop a continuum phase model, 
along the lines described by Crook et al. 179]. These authors further incorporated axonal delays 
into their model to explain the seemingly contradictory result that synchrony is stable for short 
range excitatory coupling, but unstable for long range. To see how a delay induced instability may 
arise we consider a continuum model of identical phase oscillators with space-dependent delays 

—6l(x,t) =a;-|-e J^W{x,y)H{e{y,t) - 9{x,t) - T{x,y))dy, (39) 
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where x G Z? C M and 6 G This model is naturally obtained as the continuum limit of (31) 


for a network arranged along a line with communication delays r(x, y) set by the speed of an action 
potential v that mediates the long range interaction over a distance between points in the tissue at 
X and y. Here we have used the result that for weak coupling then delays manifest themselves as 
phase-shifts. The function W sets the anatomical connectivity pattern. It is convenient to assume 
that interactions are homogeneous and translation invariant, so that W{x,y) = W{\x — y\) with 
T{x,y) = \x — y\/v, and either assume periodic boundary conditions or take H = M. 

Following m we construct travelling wave solutions of equation (|39[) for H = M of the form 


9{x,t) = Ilf + fix, with the frequency VL satisfying the dispersion relation 


— cu T c 


dyW{\y\)H{liy - \y\/v) 


When /? = 0, the solution is synchronous. To explore the stability of the travelling wave we linearise 
(39) about + fix and consider perturbations of the form to find that travelling phase 


waves solutions are stable if Re \{p) <0 for all p 7 ^ 0, where 

poo 

\{p) = e W{\y\)H'(l3y - |y|/u)[e*™ - l]dy. 


Note that the neutrally stable mode A(0) = 0 represents constant phase-shifts. For example, for 
the case that H{9) = sin 9 then we have that 


ReA(p) = 7 re[A(p,^+)-t-A(-p,/3_)], 


where 


A(p,/3) = We(p + /?) + We(p-/3)-2W,(/3), Wc{p)= W(y) cos(py)dy, 

Jo 

and fi± = ±/3 — l/(u). A plot of the region of wave stability for the choice W{y) = exp(—|y |)/2 
e > 0 in the (/3, v) plane is shown in Fig. 


18 Note that the synchronous solution /3 = 0 


and 

is unstable for small values of v. For a discussion of more realistic forms of phase interaction, 
describing synaptic interactions, see m- 


6.3.1 Phase turbulence 

For appropriate choice of the anatomical kernel and the phase interaction function, continuum 


models of the form (39) can also support weak turbulent solutions reminiscent of those seen in 


the Kuramoto-Sivashkinsky (KS) equation. The KS equation generically describes the dynamics 
near long wavelength primary instabilities in the presence of appropriate (translational, parity and 
Galilean) symmetries, and is of the form 


9t = -a9xx + fi{9xf‘ - l9x 


(40) 


where a, fi,'y > 0. For a further discussion of this model see [255]. For the model (39) with decaying 
excitatory coupling excitation and purely sinusoidal phase coupling, simulations on a large domain 
show a marked tendency to generate phase-slips and spatiotemporal pattern shedding, resulting in a 
loss of spatial continuity of 9{x, t). However, Battogtokh [38] has shown that a mixture of excitation 
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Figure 18: Stability region (black) for a phase wave 9{x, t) = hit + fix in the (/3, v) plane for the 
choice H{d) = sin 6 <, W{y) = exp(—|y|)/2 and e > 0. 


and inhibition with higher harmonics in the phase interaction can counteract this tendency and 
allow the formation of turbulent states. To see how this can arise consider an extension of (39) to 
allow for a mixing of spatial scales and nonlinearities in the form 




M 


(9(x,t) = w + / W^{x - y)H^{e{y,t) - e{x,t))dy, 

M=i 


(41) 


where we drop the consideration of axonal delays. Using the analysis of § 6.3 the synchronous wave 
solution will be stable provided \{p) <0 for all p 7 ^ 0 where 


M 


= [tu^(p)-ty^(0) 

^l=l 


WM = / 


After introducing the complex function z{x, t) 
functions as Fourier series of the form H^{9) = 


= exp(i 0 (x,t)) and writing the phase interaction 
then we can re-write (41) as 


zt = iz 


M 


UJ 


+ EE 

/i=l n^'L 


e^HUz- 




(42) 


where 

i^n{x,t)= [ W^{x-y)z"^{y,t)dy = Wf,^z'^. 

Jm 

The form above is useful for computational purposes, since ipn can easily be computed using a 
fast Fourier transform (exploiting its convolution structure). Battogtokh |38) considered the choice 
M = 3 with Hi{9) = H2{9) = sin(0 -|- a), = sin(20 -|- a) and lU^(x) = 7 ^ exp(— 7 ^|x |)/2 with 
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72 = 73 - In this case Wfj,{p) = 7 ^/( 7 ^ so that 


Kp) 


—p^ cos(a) 


f I ^2 + 263 \ 

V7?+P^ 7i+PV 


By choosing a mixture of positive and negative coupling strengths the spectrum can easily show 
a band of unstable wave-numbers from zero up to some maximum as shown in Fig. [T^ Indeed 



Figure 19; Spectrum for the phase oscillator continuum model (41) with a mixture of spatial 
scales and nonlinearities. Here Hi{6) = H2{0) = sin(0 -|- a), = sin(20 -|- a) and W^{x) = 

7 ^exp(- 7 ^|x |)/2 with 72 = 73 . Parameters are ei = 0.5, 62 = 0.15, €3 = -0.3, 71 = 1/2, 72 = 1/4, 
and a = —1.45. There is a band of unstable wave-numbers with p G (0,pc)) with pc — 1.25. 


this shape of spectrum is guaranteed when > 0 and ^M-I^u(0)/7n < 0. Similarly 


the KS equation (40) has a band of unstable wave-numbers between zero and one (with the most 
unstable wave-number at l/\/2). For the case that all the spatial scales 7 “^ are small compared 
to the system size then we may develop a long wavelength argument to develop local models for 
■f/n- To explain this we first construct the Fourier transform ^n{p,t) = W^{p)fn{p,t), where /„ is 
the Fourier transform of and use the expansion Wfj_{p) ~ 1 — 7 /l^p^ -|- 7 /l^p^ -|- .... After inverse 
Fourier transforming we find 






Noting that = H 2 = = exp(za)/(2i) = F with all other Fourier coefficients zero then |42^ 

yields 


Ot = f2-h2ReF dxxxx T ■ • •) G 

11 = 1,2 

+ €36 (73 ‘^dxx — 73 ^dxxxx + • • ■) e^*^, ( 43 ) 


where = w -|- ^ 


,,-2 


Expanding (43) to second order gives 9t = ^ — aOxx + I3{0x)^, 

Going to higher order yields fourth 

given by 7 = 


and /3 = - Cf, 




where a = l^l 

order terms and we recover an equation of KS type with the coefficient of —9x 

e^PI/( 0 ) 7 “^. To generate phase turbulence we thus require a > 0, which is also a condition 
required to generate a band of unstable wave-numbers, and /3 ,7 > 0. A direct simulation of the 
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model with the parameters for Fig. 19, for which a,l3,j > 0, shows the development of a phase 
turbulent state. This is represented in Fig. 20 where we plot the absolute value of the complex 
function 'I' = {eiWi + e 2 VF 2 ) ® z + e^Ws 0 z^. 



Figure 20: The emergence of a turbulent phase state in a phase oscillator continuum model. The 
parameters are those as in Fig. 19 with w = 0 for which a = 0.63, /3 = 5.16 and 7 = 0.096. The 
physical domain size is 2^ and we have used a numerical mesh of 2^^ points with Matlab ode45 
to evolve equation (42) with convolutions computed using fast Fourier transforms. As an order 
parameter describing the system we have chosen |T|, where T = (eiVFi + e 2 W 2 ) 'Si z + S z^. 


7 Heteroclinic attractors 


In addition to dynamically simple periodic attractors with varying degrees of clustering, the emer¬ 


gent dynamics of coupled systems such as (31) can be remarkably complex even in the case of global 


coupling. In this case, the dynamical complexity depends only on the phase interaction function 
H and the number of oscillators N. Chaotic dynamics [H] can appear in four or more globally 
coupled phase oscillators for phase interaction functions of sufficient complexity. We focus now 
on attractors that are robust and apparently prevalent in many such systems: robust heteroclinic 
attractors. 

In a neuroscience context such attractors have been investigated under several related names, 
including slow switching |135( I165| I166( 1160] where the system evolves towards an attractor that 
displays slow switching between cluster states where the switching is on a timescale determined 
by the noise, heteroclinic networks P EH [ 22 ] or winnerless competition [262] 1265] where there 
are a number of possible patterns of activity that compete with each other but such that each 
pattern is unstable to some perturbations that take it to another pattern - this can be contrasted 
to winner-takes-all competition where there is attraction to a asymptotically stable pattern. 
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These attractors obtain their dynamical structure from the presence of invariant subspaces for 
the dynamics that allow the possibility of robust connections between saddles. These subspaces 
may be symmetry-induced fixed point subspaces, spaces with a given clustering, subspaces forced 
by multiplicative coupling or subspaces forced by assumptions on the nature of the coupling. In all 
cases, there will be a number of dynamically simple nodes, usually equilibria or periodic orbits, say 
Xi, i = 1,... ,k each of which has saddle type. These nodes have unstable manifolds that, within 
some invariant subspace, limit to other nodes within the network - usually because there is a robust 
(transverse) saddle-to-sink connection between the nodes within some invariant subspace; see [29]. 
It can be verified that such heteroclinic networks can be attracting if, in some sense, the rate of 
expansion away from the nodes is not as strong as the rate of contraction towards the nodes - see 
m for some precise results in this direction. Figure illustrates some of the ingredients for a 
robust heteroclinic attractor. 



Figure 21: Schematic diagram showing a trajectory x{t) (solid line) approaching part of a robust 
heteroclinic network (bold dashed lines). The nodes Xi represent equilibria or periodic orbits of 
saddle type and the invariant subspaces Pi are forced to exist by model assumptions and there are 
connecting (heteroclinic) orbits Cj that connect the nodes within the Pi in a robust manner. A 
neighbourhood of the connecting orbits Ci is an absorbing stable heteroclinic channel that can be 
used to describe various aspects of neural system function in systems with this dynamics; see for 
example |264| . 


7.1 Robust heteroclinic attractors for phase oscillator networks 


Hansel et al. [ns] considered the dynamics of ( |31| ) with global coupling and phase interaction 
function of the form (38) for (r, a) fixed parameters. For large N, they find an open region in 
parameter space where typical attractors are heteroclinic cycles that show slow switching between 
states where the clustering is into two clusters of macroscopic size. This dynamics is examined in 
more depth in |165j where the simulations for typical initial conditions show a long and intermittent 
transient to a two-cluster state that, surprisingly, is unstable. This is a paradox because only low 
dimensional subset of initial conditions (the stable manifold) should converge to a saddle. The 
resolution of this paradox is a numerical effect: as the dynamics approaches the heteroclinic cycle 
where the connection is in a clustered subspace, there can be numerical rounding into the subspace. 
For weak perturbation of the system by additive noise, they find that the heteroclinic cycle is 
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approximated by a roughly periodic transition around the cycle whose approximate period scales 
as the logarithm of the noise amplitude. 

The bifurcations that give rise to heteroclinic attractors in this system on varying (r, a) is quite 
complex even for small N. As discussed in [27] one can only find attracting robust heteroclinic 


attractors in (31), (38) for > 4: in this case Fig. 16 shows a region where robust heteroclinic 
attractors of the type illustrated in Fig. 22 appear. A trajectory approaching such a network will 
spend much of its time near a cluster state with two groups of two oscillators each. Each time there 
is a connection between the states, one of the groups will break clustering for a short time, and over 
a longer period the clusters will alternate between breaking up and keeping together. Qualitatively 
similar heteroclinic attractors can be found for example in coupled Hodgkin-Huxley type limit cycle 
oscillators with delayed synaptic coupling as detailed in [29] and illustrated in Fig. 23 



Figure 22: A robust heteroclinic cycle attractor for the all-to-all coupled 4-oscillator system (31) 
with phase interaction function (38) and an open region of parameter space, as in m- The figure 
shows the cycle in the 3 dimensional space of phase differences; the vertices R all represent the 
fully symmetric (in-phase) oscillations (ip, ip, ip, ip), varying by 27r in one of the arguments. The 
point Q represents solutions (p, p, p n, p + tt) with symmetry (52 x 52 ) Z 2 . The heteroclinic 

cycle links two saddle equilibria Pi = (p, p, p + a, p + a) and P 2 = (p,p,p + 2 tt — a,p + 27r — a) 
with S 2 X S 2 isotropy. The robust connections Gi and G 2 shown in red lie within two dimensional 
invariant subspaces with isotropy S 2 while the equilibria 5 have isotropy S 3 . This structure is an 
attractor for parameters in the region indicated in Figure [TG] 
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Figure 23: Robust heteroclinic cycle attractor for all-to-all coupled oscillatory Hodgkin-Huxley 
systems with delay synaptic coupling. The top panel shows synchronisation indices rij that only 
equal one when the systems i and j are synchronised while the bottom panel shows the action 
potentials Vi for the four oscillators; see [29] for more details. Observe that the mechanism of 


symmetry breaking and re-synchronisation of pairs is the same as in Fig. 22 
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The heteroclinic attractors that appear for N > 4 can have more complex structures. For N = 5 
this is investigated in [241125| for (31), (|38|) and in [2Tj for more general phase interaction function 


H{ip) = — sin((/7 + a) + r sin(2(^ + /3), 


(44) 


where a, (3 and r are parameters. Figureillustrates a trajectory of (31) with global coupling and 
phase interaction function (44) as a raster plot for five phase oscillators with parameters r = 0.2, 
a = 1.8, /3 = —2.0 and oj = 1.0 and with addition of noise of amplitude 10“^. Observe there is a 
dynamic change of clustering over the course of the time-series with a number of switches taking 
place between cluster states of the type 

(01, • • • , 05) = 1 , 1, 1 , 1 ) + (y, 0,0, ^ b) 

where y, g and h represent different relative phases that are coloured “yellow”, “green” or “blue” 


in Fig. 25 to other symmetrically related cluster states. One can check that the group orbit of 
states with this clustering gives 30 symmetrically related cluster states for the system; details of 
the connections are shown in Fig. All cluster states connect together to form a single large 
heteroclinic network that is an attractor for typical initial conditions |21) . Fig. 24 illustrates the 


possible switchings between phase differences near one particular cluster state for the heteroclinic 
network in this case. The dynamics of this system can be used for encoding a variety of inputs, 
as discussed in |355j where it is shown that a constant heterogeneity of the natural frequencies 
between oscillators in this system leads to a spatio-temporal encoding of the heterogeneities. The 
sequences produced by the system can be used by a similar system with adaptive dynamics to learn 
a spatio-temporal encoded state 1212]. 

Robust heteroclinic attractors also appear in a range of coupled phase oscillator models where 
the coupling is not global (all-to-all) but such that it still preserves enough invariant subspaces for 
the connections to remain robust. For example, |154| study the dynamics of a network “motif” 
of four coupled phase oscillators and find heteroclinic attractors that are “ratchets”, i.e. they are 
robust heteroclinic networks that wind preferentially around the phase space in one direction - this 
means that under the influence of small perturbations, phase slips in only one direction can appear. 



Figure 24; Raster plot showing a robust heteroclinic attractor in a system of five globally coupled 
phase oscillators (31) with phase interaction function (44) and a particular choice of parameters 
(see text). The vertical dark lines mark times [t in horizontal axis) when the oscillator represented 
by 0fc(t) {k in vertical axis) passes through zero. Observe that most of the time there are three 
clusters. Occasionally the clustering splits and reforms different three clusters reforms, at times 
indicated by the grey bars, approximately every 70 time units. [Adapted from [2T] ] 
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Figure 25: (a) Graph of heteroclinic connections between three clnster states for a robust hetero¬ 

clinic attractor in a system of N=5 globally conpled phase oscillators. Phase interaction function 
and parameters as in Fig. 24; see m for more details of the dynamics. Each vertex represents a 
particular saddle periodic cluster state that is a permutation of the states shown in (b-d). Note 
that there are precisely two incoming and two outgoing connections at each vertex, (b-d) show the 
relative phases of the hve oscillators, indicated by the angles of the “pendula” at the vertices of a 
regular pentagon, for a sequence of three consecutive three saddle cluster states visited on a longer 
trajectory that randomly explores graph (a) in the presence of noise. Inequivalent clusters are 
characterised by different coloured “pendula” and nnmbers where yellow corresponds to 1, green to 
2 and blue to 3. The yellow cluster is stable to cluster-splitting perturbations while the bine cluster 
is unstable to such perturbations - observe that after the connection the yellow cluster becomes the 
blue cluster while the blue cluster splits up in one of two possible way. 
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7.2 Winnerless competition and other types of heteroclinic attractor 


Heteroclinic connections play an key role in describing winnerless competition dynamics |262l I265j . 
This dynamics is usually associated with firing rate model of coupled neurons of “Lotka-Volterra” 
form 

N 

Ai 
i=i 


d 

— Xi= Xi 

dt 


MjXj 


(45) 


where Xi{t) is an instantaneous firing rate of a neuron or group of neurons, Aj is a self-excitation 
term and Aij represents coupling between the Xi, though it is also possible to consider generalisa¬ 
tions |261j . These models were originally developed for predator-prey interaction of species in an 
ecosystem. Having found wide use in population ecology, they have more recently be applied to 
evolutionary game dynamics |139j . including applications in economics. Since the seminal paper of 


May |214j it has been recognised that (45) can have robust heteroclinic attractors for an open set 
of parameter choices Xi and Aij, as long as at least three species are involved. Indeed, the system 
can show a wide range of dynamics such as “winner-takes-all” equilibria where there is a stable 
equilibrium with Xi > 0 and Xj = 0 for j ^ i (i is the “winner”), “cooperative” equilibria where 
several of the Xi are non-zero as well as non-trivial periodic or chaotic dynamics. The particular 
case of “winnerless competition” gives attractors consist of chains or networks of saddles joined by 
connections that are robust because the absence of a species (or in our case, the lack of firing of 


one neuron or group of neurons) Xj = 0 is preserved by the dynamics of the model (45). 
The simplest case of this appears in a ring of iV = 3 coupled neurons with dynamics 


—xi = xi(l -xi- ax2 - Pxs), 

^X2 = X2(l - X2- ax3 - /3xi), 

^3:3 = X3(l - X3- axi - 13 x 2 ), 

for a -)- /3 > 2 and 0 < a < 1 nn, corresponding to cyclic inhibition of the neurons in one 
direction around the ring and cyclic excitation in the other direction. This “rock-paper-scissors” 
type of dynamics leads to winnerless competition that has been applied in a variety of more complex 
models of neural systems. The local behaviour near connections in such heteroclinic attractors has 
been called “stable heteroclinic channels” [264| and used to model a variety of low-level and high- 
level neural behaviours including random sequence generation, information dynamics, encoding of 
odours and working memory. We refer to the reviews [26311264( 1261] . 

Analogous behaviour has been found in a range of other coupled systems, for example [9] or 
delayed pulse-coupled oscillators [166] 1322] |32l 1159] [53l 1232] . Recent work has also considered an 
explicit constructive approach to heteroclinic networks to realise arbitrary directed networks as a 
heteroclinic attractor of a coupled cell system [22] . 

More complex but related dynamical behaviour has been studied under the names of “chaotic 
itinerancy” (see for example [152] ), “cycling chaos” [SS] > “networks of Milnor attractors” [EH] and 
“heteroclinic cycles for chaotic oscillators” [178] . It has been is suggested that these and similar 
models are useful for modelling of the functional behaviour of neural systems [328] . 

Because heteroclinic attractors are quite singular in their dynamical behaviour (averages of 
observables need not converge, there is a great deal of sensitivity of the long-term dynamics to noise 
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and system heterogeneity), it is important to consider the effect of noise and/or heterogeneities in 
the dynamics. This leads to a finite average transition time between states determined by the level 
of noise and/or heterogeneity (which may be dne to inputs to the system) and the local dynamics 
- see for example [310] . Another useful feature of heteroclinic attractors is that they allow one to 
model “input-output” response of the system to a variety of inputs. 

8 Stochastic oscillator models 

Noise is well known to play a constructive role in the neural encoding of natural stimuli, leading 
to increased reliability or regularity of neuronal firing in single neurons |206L I313j and across pop¬ 
ulations [TOO]. From a mathematical perspective it is natural to consider how noise may affect 
the reduction to a phase-oscillator description. Naively one may simply consider the addition of 
noise to a deterministic phase-oscillator model to generate a stochastic differential equation. Indeed 
models of this type have been studied extensively at the network level to understand noise-induced 
first- and second-order phase transitions, and new phenomenon such as noise-induced synchrony 
|1161122811181j or asynchrony |155| , and noise-induced turbulence |156j . We refer the reader to the 
review by Lindner |198j for a comprehensive discussion. More recently Schwabedal and Pikovsky 
have extended the foundations of deterministic phase descriptions to irregular, noisy oscillators 
(based on the constancy of the mean first return times) |279j . Ly and Ermentrout |202j and Nakao 
et al. |229j have built analytical techniques for studying weak noise forcing, and Moehlis has 
developed techniques to understand the effect of white noise on the period of an oscillator |223j . 

At the network level (global coupling) a classic paper examining the role of external noise in IF 
populations, using a phase description, is that of Kuramoto |174j . who analysed the onset of col¬ 
lective oscillations. Without recourse to a phase-reduction it is well to mention that Medvedev has 
been pioneering a phase-amplitude approach to studying the effects of noise on the synchronisation 
of coupled stochastic limit cycle oscillators [TO3im7| . and that Newhall et al. have developed a 
Fokker-Planck approach to understanding cascade-induced synchrony in stochastically driven IF 
networks with pulsatile coupling and Poisson spike-train external driven |233j . More recent work on 
pairwise synchrony in network of heterogeneous coupled noisy phase oscillators receiving correlated 
and independent noise can be found in |201j . However, note that even in the absence of synap¬ 
tic coupling, two or more neural oscillators may become synchronised by virtue of the statistical 
correlations in their noisy input streams umM- 

8.1 Phase reduction of a planar system with state-dependent Ganssian white 
noise 

For clarity of exposition let us consider the phase reduction of a planar system described by dx/dt = 
F{x)+p{x)^{t), where ^{t) is white Gaussian noise such that (C(t)) = 0 and (^(t)^(s)) = 2D6(t — s), 
where (•) denotes averaging over the realisation of and D > 0 scales the noise intensity. We employ 
an Stratonovich interpretation of the noise (such that the chain rule of ordinary calculus holds). In 
the absence of noise we shall assume that the system has a stable T-periodic limit-cycle solution, 
with a phase that satisfies d'd/dt = 1. For weak noise perturbations the state point will leave the 
cycle, though will stay in some neighbourhood, which we denote by U. To extend the notion of a 
phase off the cycle we let be a smooth function of x such that {Vx'd, F{x)) = 1 holds for any point 
X G Z7. We shall denote the other coordinate in by p, and assume that there exists a smooth 
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coordinate transformation x —)> {'d,p)- For a noise perturbed oscillator, in the new coordinates we 
have 


= 1 + hi-d, p)C{t), Ap = p) + p)^{t), 


(46) 


where we have introduced h = (Va;i?,p), / = {VxP,F) and g = {VxPiP)- Note that the full 
coordinate transformation {'&,p) = ('i9(x), p(x)) is not prescribed here, though it is common to use 
one such that p can be interpreted as some distance from cycle. Thus although results may be 
formally developed for the system (46) they cannot be directly interpreted in terms of a given 
model until the full coordinate transformation taking one from x G to (i?, p) is given. 

One can transform (46) into a stochastic phase-amplitude equation in the ltd sense, where it 
reads 


At? = 1 + p)h{^, p) + hpi-d, p)g{'d, r)] + h{'&, p)i{t), 

P) + D[g{^{d, p)h{'&, p) + gpi'd, p)g{'d, r)] + gi-ff, p)C{t), 

where the subscripts •& and p denote partial derivatives with respect to ^ and p respectively. 
Using the ltd form we may construct a Fokker-Planck equation for the time-dependent probability 
distribution Q(t,'&,p) as 


A 

m 


Q 


3 3‘^\h‘^(D^ 

[{1 + Dih,h + hpg)}Q] + 

^ [{/ + Dig,h + gpg)}Q] + + D 


d^[g^Q] 

dp^ ’ 


(47) 


with periodic boundary condition (5(t,0,p) = Q{t,2'K, p). 

When D = 0 the steady state distribution is given by Qo{'d,p) = 6{p)/{2tt). For small noise 
amplitude D we expect Qq to still localise near p = 0 |359j . say over a region —pc < p < Pc- In this 
case it is natural to make the approximation that for large t, Q = 0 = dQjdp ai p = Fpc- Now 
introduce the marginal distribution 


P{t,'d)= / Q(t,'(?,p)dp. 


'-Pc 


Following |359j we can integrate (47) over the interval I = [—PoPc] and generate a Fokker-Planck 
equation for P. The last three terms in (47) vanish after integration due to the boundary conditions, 
so that we are left with 

r\ Cl /* d2 /* 

j^^(h^^ + K9)Q<ip + D-^ j^h'^Qdp. 


We now expand about /? = 0 to give h{'d,p) = Z{'&) + php{'d,d) + ... and g{'d,p) = g{'d,0) + 
pgp{'d,0) where Z{i!}) is identified as the infinitesimal phase response {'^x'&,p\p=o)■ In the 

limit of small D where Qo ~ 5{p)/{2'k) we note that, for an arbitrary function R{'0), 
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Making use of the above gives the Fokker-Planck equation as 

where Y{'&) = hp{'i},0)g{'d,0). Hence, the corresponding Ito equation is 

^/ = l + D[Z{i})Z'{i!}) + y(^)] + Z{mt), 
while the Stratonovich version is 


(48) 


(49) 


—^ = l + DY{'d) + Z{'d)^{t). 


(50) 


Equations (49) and (50) are the stochastic phase oscillator descriptions for a limit cycle driven by 
weak white noise. These make it clear that naively adding noise to the phase description misses 
not only a multiplication by the iPRC but also the addition of a further term Y ('d) that contains 
information about the amplitude response of the underlying limit cycle oscillator. 

We are now in a position to calculate the steady state probability distribution Po{'&) and use 
this to calculate the moments of the phase dynamics. Consider (48) with the boundary condition 
P{t,0) = P{t,2TT), and set Pt = 0. Adopting a Fourier representation for Pq, Z and Y as Po{'d) = 
Yhn , Z{'d) = Yin , T('d) = Yn^nP'^^ , allows US to obtain a set of equations for the 

unknown amplitudes Pi as 

-Pi + D ^ ZnZmi{l-n)Pi_(^n+rn)-D'^YnPi-n = Kdifl, / E Z, (51) 

for some constant K. For D = 0 we have that Pq = K, and after enforcing normalisation we may 
set K = l/(27r). For small D we may then substitute Pq into (51) and work to next order in D to 
obtain an approximation for the remaining amplitudes, / 7 ^ 0 , in the form 




ZnZmim - Yi 

{n,m |n+m=/} 


Using this we may reconstruct the distribution P{'d), for small D, from (18) as 


pm = ^ ^ imz'm - Y{tt )+ To) 


1 


c 2 tt 


Y = Y{'d)d'd. 

--1 d, 


(52) 


The mean frequency of the oscillator is defined by w = lim'r^oo P ^ Jq a 7 'f^(^)dt- This can be 


or an ar 


aitrary 27r-periodic 


calculated by replacing the time average with the ensemble average. 

function R we set limr_^oo T~^ -R('f^)-Po(i^)d'd. Using (52) and (49) we obtain 

u = l + DYo + 0{D), 

where we have used the fact that {Z{'d)^{t)) = {Z {-&)) {^{t)) = 0. We may also calculate the 
phase-diffusion D as 


D = 


f 

D 

TT 


^'d{t + t) — ( -^1? 
dt ^ ^ \dt 


dt \ dt 


dr 


f27r 


Z^{^)d^ + OiD"^), 


where we use the fact that (^(t)) = 0 and {C{t)^{s)) = 2D6{t — s). This recovers a well known 
result of Kuramoto m- 
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8.2 Phase reduction for noise with temporal correlation 


A recent paper by Teramae et al. [318] shows that when one considers noise described by an 
Ornstein-Uhlenbeck (OU) process with a finite correlation time then this can interact with the 
attraction time-scale of the limit cycle and give fundamentally different results when compared to 
Gaussian white noise (which has a zero correlation time). This observation has also been indepen¬ 
dently made in |358| . Both approaches assume weak noise, though |358) makes no assumptions 
about relative time-scales, and is thus a slightly more general approach than that of |318j . Re¬ 
lated work by Goldobin et al. for noise r]{t) with zero-mean {r]{t)) = 0 and prescribed 

auto-correlation function C'(r) = {r]{T)r]{0)) , yields the reduced Stratonovich phase equation 


where 


—^ = l + DY{^) + Z{^)r]{t), 
at 


~ 1 

Yi^) = 0) g{t} - V’, 0)C(V')e-^^dV^, 


(53) 


(54) 


where A is the average rate of attraction to the limit cycle. Note that for C(t) = 2D5{t), Y{'d) = 
y('d) and (53) reduces to (50) as expected. To lowest order in the noise strength the steady state 
probability distribution will simply be Po{'&) = l/(27r). Therefore to lowest noise order the mean 
frequency is determined from an ensemble average as 


w = w + DYo + 


1 

47r 


(*2n PC 

D Jo 


dV’^(??-V’)G(V'), 


where the last term comes from using the Ito form of (53) and the subscript 0 notation is defined 
as in (52). The phase-diffusion coefficient at lowest noise order is given by 

POO 

di? / drZ(i?)Z(79 + r)C(r). 


D = — 
2 -k 


r-27r 


Let us now consider the example of OU noise so that C'(r) = D'j exp{—'y\T\). Furthermore 
let us take the simultaneous limit 7 —>■ oo (zero correlation timescale) and A —)• 00 (infinitely fast 


attraction), such that the ratio a = A/7 is constant. In this case we have from (54) that 


Y{^) = 


Yj'd) 
\ (y. 


(55) 


Hence, when the correlation time of the noise is much smaller than the decay time constant a = 0 
and we recover the result for white Gaussian noise. In the other extreme when a —>■ 00, where 
the amplitude of the limit cycle decays much faster than the correlation time of the noise, then Y 
vanishes and the reduced phase equation is simply d'd/dt = 1 -|- Z{'&)r](t), as would be obtained 
using the standard phase reduction technique without paying attention to the stochastic nature of 
the perturbation. 


9 Low dimensional macroscropic dynamics and chimera states 

The self-organisation of large networks of coupled neurons into macroscopic coherent states, such 
as observed in phase-locking, has inspired a search for equivalent low-dimensional dynamical de¬ 
scriptions. However, the mathematical step from microscopic to macroscopic dynamics has proved 
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elusive in all but a few special cases. For example, neural mass models of the type described in 


2.6 only track mean activity levels and not the higher order correlations of an underlying spiking 


model. Only in the thermodynamic limit of a large number of neurons firing asynchronously (pro¬ 
ducing null correlations) are such rate models expected to provide a reduction of the microscopic 
dynamics. Moreover, even here the link from spike to rate is often phenomenological rather than 
rigorous. Unfortunately only in some rare instances has it been possible to analyse spiking networks 
directly (usually under some restrictive assumption such as global coupling) as in the spike-density 
approach [238], which makes heavy use of the numerical solution of coupled PDFs. Recently how¬ 
ever, exact results for globally pulse-coupled oscillators described by the Winfree model |352| have 
been obtained by Pazo and Mont brio |249j . This makes use of the Ott-Antonsen (OA) ansatz, 
which was originally used to find solutions on a reduced invariant manifold of the Kuramoto model 
m- The major difference between the two phase-oscillator models being that the former has 
interactions described by a phase product structure and the latter a phase difference structure. 


9.1 Ott-Antonsen reduction for the Winfree model 


The Winfree model is described in § [^ as a model for weakly globally pulse-coupled biological 
oscillators, and can support incoherence, frequency locking, and oscillator death when P{9) = 
1 -|- cos 9 and R{9) = — sin0 (TT]. We note however that the same model is exact when describing 
nonlinear IF models described by a single voltage equation, and that in this case we do not have 
to restrict attention to weak-coupling. Indeed the OA ansatz has proved equally successful in 
describing both the Winfree network with a sinusoidal PRC and a network of QIF neurons 


[200] . This is perhaps not surprising since the PRC of a QIF neutron can be computed using (26), 
and for the case described by ([^ with r = 1 and f{v) = and infinite threshold and reset then 
R{9) = a{l — cos9) with a = I/'/i for 9 G [0, 27r) (which is the shape expected for an oscillator 
near a SNIC bifurcation). We shall now focus on this choice of PRC and a pulsatile coupling that 
we write in the form 


E' 

mez 


imd 


P{9) = 27r 5{9 — Inn) = 

nGZ 

where we have introduced a convenient Fourier representation for the periodic function P. We 
now consider the large N limit in (30) and let p{9\u},t)d9 be the fraction of oscillators with phases 
between 9 and 9 + d9 and natural frequency u at time t. The dynamics of the density p is governed 
by the continuity equation (expressing the conservation of oscillators) 


-~ v{9,t) — uj + eah{t)[l — + c *^)/2], (56) 

where the mean-field drive is h{t) = \\m.p^^^'^-N~^P{9j). Boundary conditions are periodic 
in the probability flux, namely /9(0|a;, t)u(0, t) = \\m.o^ 2 iT p{9\'^-,t)v{9,t). A further reduction in 
dimensionality is obtained for the choice that the distribution of frequencies is described by a 
Lorentzian g{ijo) with 

7r(a;-u;o)2 +A2’ 

for fixed A and cuq (controlling the width and mean of the distribution respectively), which has 
simple poles at uj± = ujq ± iA. 
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A generalised set of order parameters is defined as 




/ oo /* 27 r 

da;5f(a;) / (16p{e\uj,t) 

■OO J 0 


AmO 


m G N. 


The integration over oj can be done using a contour in the lower half complex plane so that Zm (t) = 

p{9\ijj-,t)), where we have introduced the inner product (n(6<), v{6)) = d6u{6)v{6)d9. 

The OA ansatz assumes that the density p can be written in a restricted Fourier representation as 


27rp{9\uj, t) = 1 + ^ a{uj, + cc, 


(57) 


m=l 


where cc stands for complex conjugate. Substitution into the continuity equation (56) and balancing 
terms in e*™® shows that a must obey 


d h 

—a = + eah)a + iea—(1 + a^). 

cji 2 


(58) 


Moreover, using the inner product structure of Zm we easily see that Zm{t) = [a(a;_,t)]™. Thus 
the Kuramoto order parameter Zi = Re~^'^ is governed by (58) with w = uo- yielding: 


—i? = —AR — ea—(1 — R^) sin 
dt 2 ^ ^ ’ 


— 'I' = Wo + eah 
dt 


1 1 + ^^ t' 

1-cos 'h 


2R 


(59) 


To calculate the mean-field drive h we note that it can be written as 

r2TT 


h = 


'0 


d9p{9\uj_,t)P{9) = ^Z„ 


1 + {z,r + = 1 + ytV + i3T’ 

m=l ^ 


Hence, we have explicitly that h = h{R, T) with 

1 _ d 2 

h{R,d!) = 


1 — 2R cos T + ’ 


d<R<l. 


(60) 


The planar system of equations defined by (59) and (60) can be readily analysed using numerical 
bifurcation analysis. 

We note that the OA density (57) can be written in the succinct form 

1 -|- ae*® 


2-Kp{9\ujR) = ^ 


1 — ae 


iO 


+ cc 


The form of this equation suggests recasting the density in terms of new variables VF = (1 + 
a)/(l — a) and V = i{l + e*^)/(l — e*®) = — cot(0/2) |275j . In this case we find that p(H|w, t) has 
a Lorentzian shape given explicitly by 


1 


p{V\uj,t) = - 


u 


TT {V — vY + V? ' 


W = n + in. 
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with dV/dt = I + eh, where we have set oj = 2 ^/I and re-scaled V hy a factor of \/7- Using 
the continuity equation for p[V\oj,t) and then integrating over V gives the evolution of W as 

BW 

— =i[/ + e/r-lU2], v = j dVpiV\uj,t)V, (61) 


where we identify v as the average of V. The firing probabilities for arbitrary cells at time t 
are equal to the passage rate of the probability density through the spike phase, which gives the 
population firing rate as 

= lim p{V\u;,t)^V = -u{uj,t). 

V^oo dt TT 


Thus we may use (61) to determine the evolution of the coupled rate and average network activ¬ 
ity by integrating over the distribution of parameters I = uP’jd. The evolution of {r(t),v{t)) = 
f dujg(uj)(r(uj, t),v(uj, t)) = (r(a;o — iA, t), u(a;o — iA, t)) is then found as 


d wqA 
r = —:-h 2rv, 


dt 


2tt 


d 2 

v = v^ + - 

dt 4 


2 2 

-|- evrr — nr, 


(62) 


where we have used the result that h = {W + lU)/2 = nr. Exploiting the fact that the order 
parameters in the ‘0’ and ‘U’ descriptions of macroscopic dynamics are related by a conformal 
mapping, namely Zi = (W — 1)/(1U -|- 1), we may explore both the firing rate and degree of 
synchrony of the network. Thus ( |62[ ) is a mean field reduction of a population of spiking QIF 
neurons, and unlike a phenomenological neural mass model it is able to track a measure of within 


population synchrony. The system (62) is capable of supporting bistable fixed point behaviour as 
shown in Fig. 26, where we also plot the basin boundary (stable manifold of a saddle). The OA 



Figure 26: A phase portrait for the macroscopic ‘U’ dynamics showing the rate and average (r, u) 
plane with three fixed points for wq = 1 and A = 10 = e. For this parameter set the system 
is bistable and the stable (blue) and unstable (green) manifolds of the saddle are plotted. The 
remaining curves are nullclines. 


ansatz has also proved remarkably useful in understanding non-trivial solutions such as chimera 
states (where a sub-population of oscillators synchronises in an otherwise incoherent sea). 
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9.2 Chimera states 


Phase or cluster synchronised states in systems of identical coupled oscillators have distinct limita¬ 
tions as descriptions of neural systems where not just phase but also frequency clearly play a part 
in the processing, computation and output of information. Indeed, one might expect that for any 
coupled oscillator system that is homogeneous (in the sense that any oscillators can be essentially 
replaced by any other by a suitable permutation of the oscillators), the only possible dynamical 
states are homogeneous in the sense that the oscillators behave in either a coherent or an inco¬ 
herent way. This expectation however is not justified - there can be many dynamical states that 
cannot be easily classified as coherent or incoherent, but that seem to have a mixture of coherent 
and incoherent regions. Such states have been given the name “Chimera state” by Abrams and 
Strogatz [Sill] and have been the subject of intensive research over the past five years. For reviews 
of chimera state dynamics we refer the reader to [18311247] . 

Kuramoto and Battogtokh |176L I175j investigated continuum systems of oscillators of the form 

= to — € [ G(x — x') sm(6(x,t) — 9{x',t) + a) dx', (63) 

ot Jd 


where 9 represent phases at locations x G C M, the kernel G{u) = k exjp{— k \ u \) /2 represents 
a non-local coupling and uj,a,K are constants. Interestingly this model is precisely in the form 
presented in § 6.3 as equation (39) for an oscillatory model of cortex, although here there are 
no space-dependent delays and the interaction function is H{9) = sin(0 — a). Kuramoto and 
Battogtokh found for a range of parameters near a = '7r/2, and carefully selected initial conditions, 
that the oscillators can split into two regions in x, one region which is frequency synchronised (or 
coherent) while the other region shows a nontrivial dependence of frequency on location x. An 
example of a Chimera state is shown in Fig. 


27 



Figure 27: A snapshot of a Chimera state for the model (63) in a system of length 4 using 2® 
numerical grid points and periodic boundary conditions. Here a; = 0, e = 0.1, a = 1.45 and k = 1. 
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Note that a discretisation of (63) to a finite set of N coupled oscillators is 


d 

—t) = w - e ^ Kij sm{9i - 9j + a) 

k=l 


(64) 


where 9i G [0, 27r) represents the phase at location i = 1 ,..., and the coupling matrix Kij = 
G{\i—j\/N)/N is the discretised interaction kernel (assuming a domain of length 1). Using different 
kernels, G{u) = exp(—Kcos(27rn)) and an approximation G{u) = 1 — kcos{2ttu) for small k, Abrams 
and Strogatz [3] identified similar behaviour and [1] discussed a limiting case of parameters such 
that the continuum system provably has chimera solutions. The OA reduction discussed in § 9.1 


allows an exact reduction of oscillator networks of this form and in the continuum limit this can give 
a solvable low-order system whose solutions include a variety of chimera states |183j . It is useful to 
note that when a. = 7r/2, pure cosine coupling results in an integrable Hamiltonian system [343| . 
such that disordered initial states will remain disordered. Thus a determines a balance between 
spontaneous order and permanent disorder. 

However, it seems that chimera states are much more “slippery” in finite oscillator systems 
than in the continuum limit. In particular, Wolfrum and Omel’chenko |354| note that for finite 


approximations of the ring (63) by oscillators, with a mixture of local and nearest rA^-neighbour 
coupling corresponding to (64) with a particular choice of coupling matrix chimera states 


apparently only exist as transients. However, the lifetime of the typical transient apparently grows 
exponentially with N. Thus, at least for some systems of the form (64), chimeras appear to be 
a type of chaotic saddle. This corresponds to the fact that the boundaries between the regions 
of coherent and incoherent oscillation fluctuate apparently randomly over a long timescale. These 
fluctuations lead to wandering of the incoherent region as well as change in size of the region. 
Eventually these fluctuations appear to result in typical collapse to either fully coherent or fully 
incoherent oscillation |354j . 


Although this appears to be the case for chimeras for (64), there are networks such as coupled 


groups of oscillators; m or two dimensional lattices [289j where chimera attractors can appear. 
It is not clear what will cause a chimera to be transient or not, or indeed exactly what types of 
chimera-like states can appear in finite oscillator networks. A suggestion of [2B] is that robust 
neutrally stable chimeras may be due to the special type of single-harmonic phase interaction 
function used in ( 63|64 ). 

More recent work includes investigations of chimeras (or chimera-like states) in chemical |323) 
or mechanical oscillator networks |210j ; chimeras in systems of coupled oscillators other than phase 
oscillators have been investigated in many papers; for example in Stuart-Landau oscillators [1761 
12411 I284j . Winfree oscillators |249j and models with inertia |239j . Other recent work includes 
discussion of feedback control to stabilise chimeras [292j . investigations of chimeras with multiple 
patches of incoherence |207j and multicluster and traveling chimera states |357j . 

In a neural context chimeras have also been found in pulse-coupled LIE networks [240) . and 
hypothesised to underly coordinated oscillations in unihemispheric slow-wave sleep, whereby one 
brain hemisphere appears to be inactive while the other remains active jM!- 


10 Applications 

We briefly review a few examples where mathematical frameworks are being applied to neural 
modelling questions. These cover functional and structural connectivity in neuroimaging, central 
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pattern generators (CPGs) and perceptual rivalry. There are many other applications we do not 
review, for example to deep brain stimulation protocols [203] or to modelling of epileptic seizures 
where network structures play a key role |319j . 


10.1 Functional and structural connectivity in neuroimaging 


Functional connectivity (FC) refers to the temporal synchronisation of neural activity in spatially 
remote areas. It is widely believed to be significant for the integrative processes in brain function. 
Anatomical or structural connectivity (SC), is widely believed to play an important role in deter¬ 
mining the observed spatial patterns of FC. However, there is clearly a role to be played by the 
dynamics of the neural tissue. Even in a globally connected network we would expect this to be the 
case, given our understanding of how synchronised solutions can lose stability for weak coupling. 
Thus it becomes useful to study models of brain like systems built from neural mass models (such 
as the Jansen-Rit model of § 2.6), and ascertain how the local behaviour of the oscillatory node 
dynamics can contribute to global patterns of activity. For simplicity, consider a network of N 
globally coupled identical Wilson-Cowan |349j neural mass models: 


dt" 


-Xi = —X-, 


+ P + Cif{xi) - C2f{yi) 


N 


-TiVi = -yi + Q + csf{xi) - Cifiyi), 

at 

for i = 1,..., N. Here {xi,yi) G represents the activity in each of a pair of coupled neuronal 
population models, / is a sigmoidal firing rate given by /(x) = 1/(1 -|- e“*) and {P,Q) represent 
external drives. The strength of connections within a local population is prescribed by the co¬ 
efficients Cl,... C 4 , which we choose as ci = C 2 = cs = 10 and C 4 = —2 as in |142j . For e = 0 it is 
straight forward to analyse the dynamics of a local node and find the bifurcation diagram in the 
(P, Q) plane as shown in Fig. 28 Moreover, for e <C 1 we may invoke weak coupling theory to 


describe the dynamics of the full network within the oscillatory regime bounded by the two Hopf 


curves shown in Fig. 28 left. From the theory of § 6.2.1 we would expect the synchronous solution 
to be stable if €p['{0) > 0. Taking e > 0 we can consider P['{0) as a proxy for the robustness 
of synchrony. The numerical construction of this quantity, as in ra, predicts that there will be 
regions in the {P,Q) plane associated with a breakdown of FC (where H'{0) < 0 ), as indicated 
by points a and b in Fig. This highlights the role that local node dynamics can have on 
emergent network dynamics. Moreover, we see that simply by tuning the local dynamics to be 
deeper within the existence region for oscillatory solutions we can, at least for this model, enhance 
the degree of FC. It would be interesting to explore this simple model further for more realistic 
brain like connectivities, along the lines described in m- Moreover, given that this would preclude 
the existence of the synchronous state by default (since we would neither have that P[{0) =0 or 
Wij would be independent of i) then it would be opportune to explore the use recent ideas in 
[23511251j to understand how the system could organise into a regime of remote synchronisation 
whereby pairs of nodes with the same network symmetry could synchronise. For related work on 
Wilson-Cowan networks with some small dynamic noise see ISD, though here the authors construct 
a phase-oscillator network by linearising around an unstable fixed point, rather than use the notion 
of phase response. 
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Figure 28: Left: Bifurcation diagram for an isolated Wilson-Cowan node in the (P, Q) plane. Here 
HB denotes Hopf bifurcation and SN a saddle node of fixed-points bifurcation. At points a and c 
we find H'(0) < 0 and at point b H'{0) > 0. A breakdown of FC (loss of global synchrony) within 
a globally coupled network is predicted at points a and c. 


10.2 Central Pattern Generators 


CPGs are (real or notional) neural subsystems that are implicated in the generation of spatio- 
temporal patterns of activity [226j , in particular for driving the relatively autonomous activities such 


as as locomotion [68l 128711278| or for driving involuntary activities such as heartbeat, respiration 
or digestion m- These systems are assumed to be behind the creation of the range of walking 
or running patterns (gaits) that appear in different animals [67]. The analysis of phase-locking 
provides a basis for understanding the behaviour of many CPGs, and for a nice overview see the 
review articles by Harder and Bucher [209] and Hooper |141j . 

In some cases, such as the Leech [Hirudo medicinalis) Heart or Caenorhabditis elegans locomo¬ 
tion, the neural circuitry is well studied. For more complex neural systems and in more general 
cases CPGs are still a powerful conceptual tool to construct notional minimal neural circuitry 
needed to undertake a simple task. In this notional sense they have extensively been investigated 
to design control circuits for actuators for robots; see for example the review [145] . Recent work 
in this area includes robots that can reproduce salamander walking and swimming patterns m- 
Since the control of motion of autonomous “legged” robots is still a very challenging problem in 
real-time control, one hope of this research is that nature’s solutions (for example, how to walk 
stably on two legs) will help inspire robotic ways of doing this. 

CPGs are called upon to produce one or more rhythmic patterns of actuation; in the particular 
problem of locomotion, a likely CPG is one that will produce the range of observed rhythms of 
muscle actuation, and ideally the observed transitions between then. For an early discussion of 
design principles for modelling CPGs, see m- This is an area of modelling where consideration of 
symmetries as in § 3.6 has been usefully applied to constrain the models. For example |69j examine 


models for generating the gaits in a range of vertebrate animals, from those with two legs (humans) 
through those with four (quadrupeds such as horses have a wide range of gaits - walk, trot, pace. 
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(a) 
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Figure 29: Schematic diagram showing central pattern generators (a) with 4n coupled cells that is 
used to model gait patterns in animals with 2n legs |124j and (b) a three cell motif of bursters with 
varying coupling strengths, as considered in |288[ 1353] . 


canter, gallop - they may use) or larger numbers of legs (myriapods such as centipedes). Insects 
make use of six legs for locomotion while other invertebrates such as centipedes and millipedes 
have a large number of legs that are to some extend independently actuated. As an example, 
|124j consider a schematic CPG of 2n oscillators for animals with n legs, as shown in Fig. [2^a). 


The authors use symmetry arguments and Theorem 3.1 to draw a number of model-independent 
conclusions from the CPG structure. 

One can also view CPGs as a window into more fundamental problems of how small groups of 
neurons coordinate to produce a range of spatio-temporal patterns. In particular, it is interesting 
to see how the observable structure of the connections influences the range and type of dynamical 
patterns that can be produced. For example, |288) consider a simple three-cell “motif” networks 
of bursters and classify a range of emergent spatio-temporal patterns in terms of the coupling 
parameters. Detailed studies |353| investigate properties such multistability and bifurcation of 
different patterns and the influence of inhomogeneities in the system. This is done by investigating 
return maps for the burst timings relative to each other. 

The approach mentioned in 


3.11 and outlined in |307l I123j seems to be a promising way of 


providing a context in which view GPGs where there is a constrained connection structure suggested 
by either the neurobiology or from conceptual arguments, but this structure is not purely related 
to symmetries of the network. For example, use that formalism to understand possible spatio- 
temporal patterns that arise in lattices or [151j that relates synchrony properties of small motif 
networks to spectral properties of the adjacency matrix. 


10.3 Perceptual rivalry 

Many neural systems process information - they need to produce outputs that depend on inputs. If 
the system effectively has no internal degrees of freedom then this will give a functional relationship 
between output and input so that any temporal variation in the output corresponds to a temporal 
variation of the input. However, this is not the case for all but the simplest systems and often 
outputs can vary temporally unrelated to the input. A particularly important and well-studied 
system that is a model for autonomous temporal output is perceptual rivalry, where conflicting 
information input to a neural system results, not in a rejection or merging of the information, but 
in an apparently random “flipping” between possible “rival” states (or percepts) of perception. 
This nontrivial temporal dynamics of the perception appears even in the absence of a temporally 
varying input. The best studied example of this type is binocular rivalry, where conflicting inputs 
are simultaneously made to each eye. It is widely reported by subjects that perception switches 
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from one eye to the other, with a frequency that depends on a number of factors such as the 
contrast of the image [H]. More general perceptual rivalry, often used in “optical illusions” such 
as ambiguous figures - the Rubin vase, the Necker cube - show similar behaviour with percepts 
shifting temporally between possible interpretations. 

Various approaches |290j have been made to construct nonlinear dynamical models of the gen¬ 
eration of a temporal shifting between possible percepts such as competition models [15], bifur¬ 
cation models, ones based on neural circuitry [182j . or conceptual ones [350] based on network 
structures [inj or on heteroclinic attractors |20j . 

11 Discussion 

As with any review we have had to leave out many topics that will be of interest to the reader. 
In particular we have conhned ourselves to “cell” and “system-level” dynamics rather that “sub- 
cellular” behaviour of neurons. We briefly mention some other active areas of mathematical research 
relevant to the science of rhythmic neural networks. Perhaps the most obvious area that we have not 
covered in any depth is that of single unit (cell or population) forcing, which itself is a rather natural 
starting point for gaining insights into network behaviour and how best to develop mathematical 
tools for understanding response |221( 1184] . For a general perspective on mode-locked responses to 
periodic forcing see |115j and |255j . 

For a more recent discussion of the importance of mode-locking in auditory neuroscience see 
|186l 1190] and in motor systems see |80j . However, it is well to note that not much is known 
about nonlinear systems with three or more interacting frequencies |62|, as opposed to periodically 
forced systems where the notions of Farey tree and the devil’s staircase have proven especially 
useful. We have also painted the notion of synchrony with a broad mathematical brush, and not 
discussed more subtle notions of envelope locking that may arise between coupled bursting neurons 
(where the within burst patterns may desynchronise) jllj . This is especially relevant to studies 
of synchronised bursting |282j and the emergence of chaotic phenomena |227j . Indeed, we have 
said very little about coupling between systems that are chaotic, such as described in |361j . the 
emergence of chaos in networks |297[ 1296] or chaos in symmetric networks [H] . 

The issue of chaos is also relevant to notions of reliability, where one is interested in the stability 
of spike trains against fluctuations. This has often been discussed in relation to stochastic oscillator 
forcing rather than those arising deterministically in a high dimensional setting |321[ 11181 llOOl 1195] . 
Of course, given the sparsity of firing in cortex means that it may not even be appropriate to treat 
neurons as oscillators. However, some of the ideas developed for oscillators can be extended to 
excitable systems, as described in [1441 1266] . As well as this it is important to point out that 
neurons are not point processors, and have an extensive dendritic tree, which can also contribute 
significantly to emergent rhythms as described in |47l I28U] . as well as couple strongly to glial cells. 
Although the latter do not fire spikes, they do show oscillations of membrane potential m- At the 
macroscopic level it is also important to acknowledge that the amplitude of different brain waves 
can also be significantly affected by neuromodulation [188] . say through release of norepinephrine, 
serotonin and acetylcholine, and the latter is also thought to be able to modulate the PRC of a 
single neuron [309] . 

This review has focused mainly on the embedding of weakly coupled oscillator theory within a 
slightly wider framework. This is useful in setting out some of the neuroscience driven challenges 
for the mathematical community in establishing inroads into a more general theory of coupled 
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oscillators. Heterogeneity is one issue that we have mainly side-stepped, and remember that the 
weakly coupled oscillator approach requires frequencies of individual oscillators to be close. This 
can have a strong effect on emergent network dynamics |335] . and it is highly likely that even a 
theory with heterogeneous phase response curves |327| will have little bearing on real networks. 
The equation-free coarse-graining approach may have merit in this regard, though is a numerically 
intensive technique [^ . 

We suggest a good project for the future is to develop a theory of strongly coupled heterogeneous 
networks based upon the phase-amplitude coordinate system described in § 5.2, with the challenge 


to develop a reduced network description in terms of a set of phase-amplitude interaction functions, 
and an emphasis on understanding the new and generic phenomena associated with nontrivial 
amplitude dynamics (such as clustered phase-amplitude chaos and multiple attractors). To achieve 
this one might further tap into recent ideas for classifying emergent dynamics based upon the group 
of structural symmetries of the network. This can be computed as the group of automorphisms 
for the graph describing the network. For many real-world networks, this can be decomposed into 
direct and wreath products of symmetric groups [205]. This would allow the use of tools from 
computational group theory |140j to be used, and open up a way to classify the generic forms of 
behaviour that a given network may exhibit using the techniques of equivariant bifurcation theory. 


Appendix A 

The Hodgkin-Huxley description of nerve tissue is completed with: 

am{V) = ^ , MV) = 0.07exp[-0.05(H + 65)], 

MV) = 1 _ expf-IuF + 55)] ’ /5m(H) =4.0exp[-0.0556(H + 65)], 

MV) = i^,,p[_o\(^ + 35)] ’ = 0-125exp[-0.0125(H + 65)], 

and C = 1//F cm“^, gi = 0.3mmho cm“^, gx = 36mmho cm“^, gxa = 120mmho cm“^, Vl = 
—54.402mV, Vk = —77mV and Vxa = 50mV. (All potentials are measured in mV, all times in ms 
and all currents in gA per cm^). 


Glossary 

We give a brief list of some of the abbreviations used within the review. 

DDE Delay differential equation 

IF Integrate and fire (model for neural oscillator) 

iPRC Infinitesimal phase response curve 

ISI Inter-spike interval 

FHN Fitzhugh-Nagumo equation (model for neural oscillator) 
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HH Hodgkin-Huxley equation (model for neural oscillator) 

LIF Leaky integrate and fire (model for neural oscillator) 

ML Morris-Lecar equation (model for neural oscillator) 

MSF Master stability function 
ODE Ordinary differential equation 
PDE Partial differential equation 
PRC Phase response curve 

QIF Quadratic integrate and fire (model for neural oscillator) 

SDE Stochastic differential equation 
SHC Stable heteroclinic channel 

SNIC Saddle-node on an invariant circle (bifurcation) 

WLC Winnerless competition 
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